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1.  INTRODUCTION  AND  SUMMARY 


This  report  summarizes  the  work  performed  on  the  "Time  Domain  Algorithms" 
project,  under  contract  No.  F30602-84-C-0016.  The  objective  of  this  research 
was  to  address  some  of  the  Issues  related  to  the  non-statlonary,  time-varying 
nature  of  signals  propagating  through  a  communications  channel. 

Many  communications  signals  (such  as  frequency  modulation  and  phase 
modulation  formats)  are  characterized  by  a  constant  modulus  property.  The 
propagation  of  these  signals  through  the  communication  channel  Introduces 
various  types  of  “distortions"  due  to  mil tl path,  fading  and  similar  effects. 
These  distortions  tend  to  destroy  the  constant  modulus  properties  of  the 
signal.  Recently,  an  adaptive  processing  technique  was  developed  which  has 
the  capability  of  correcting  the  effects  of  the  channel  and  restoring  the 
constant  modulus  property  of  the  signal.  This  technique  has  been  shown  to  be 
very  useful  In  communications  applications.  We  have  analyzed  this  Important 
processing  technique  and  developed  some  useful  extensions. 

In  Appendix  A  we  prove  global  convergence  of  the  Constant-Modulus 
Algorithm  (CMA)  for  the  case  of  a  real  channel  when  the  model  order  Is  equal 
to  or  greater  than  that  of  the  channel  (the  so-called  "model -complete" 
case).  The  analysis  Is  based  on  an  exact  fourth-order  Taylor  reries 
representatl on  for  the  cost  function  minimized  asymptotically  by  the  CMA. 

In  Appendix  B  we  present  several  extensions  to  the  CMA  Including  IIR 
equalization,  a  real-signal  version  having  properties  as  good  as  the  complex 
version,  use  of  the  Gauss-Newton  method  In  place  of  gradient  descent. 
Interference  rejection,  and  more.  Some  preliminary  simulation  results  are 
presented  In  Appendix  F. 

The  processing  and  estimation  of  signals  propagating  through  time-varying 
channels  requires  time- varying  filter  structures.  While  the  area  of  time- 
invariant  digital  filters  Is  well  developed,  relatively  little  fundamental 
work  Is  available  on  the  time-varying  case.  We  Investigated  some  of  the  basic 
questions  related  to  time-varying  digital  filtering  and  were  able  to  derive  a 
novel  filter  structure  for  such  applications. 


In  Appendix  C  we  derive  the  set  of  finite-order,  linear,  time-invariant 
filters  by  sampling  lossless  propagation  through  a  variable-impedance 
medium.  This  leads  to  a  flexible  class  of  time-varying  filter  structures, 
termed  "Waveguide  Filters''  (WGF)  In  which  signal  power  Is  decoupled  from 
changes  In  the  filter  parameters.  These  structures  are  "balanced"  In  the 
sense  that  the  decoupling  between  signal  power  and  time-varying  filter 
coefficients  is  maintained  for  each  individual  section  In  the  structure.  In 
addition,  limit  cycles  and  overflow  oscillations  are  suppressed,  even  In  the 
time-varying  case,  when  implemented  with  "passive"  arithmetic.  Finally,  the 
WGF  structures  can  be  Interconnected  In  series  or  in  parallel  in  a  way  which 
does  not  disturb  the  signal /coefficient  decoupling  or  the  power  balance. 

Thus, the  waveguide  filters  are  very  useful  for  modeling  physical  systems,  and 
the  exactness  of  their  physical  Interpretation  enhances  their  suitability  for 
the  time-varying  case.  All  results  are  obtained  for  the  multi-input/multi- 
output  case. 

Another  topic  which  was  briefly  Investigated  during  this  study  was  the 
adaptive  equalization  of  rapidly  time-varying  multipath  channels.  Due  to  time 
limitations  only  a  conceptual  study  of  this  topic  was  possible.  In  Appendix  0 
we  describe  on  adaptive  equalizer  for  eliminating  distortion  due  to  multipath 
propagation  in  high-speed  digital  radio  systems.  The  equalizer  Is  aimed  at 
the  case  of  very  fast  channel  fluctuation,  where  "fast"  Is  defined  relative  to 
the  Impulse-response  duration  of  the  Inverse  of  the  Instantaneous  multipath 
transfer  function.  The  method  Is  of  the  decision- feedback  type  where  the 
demodulated  symbols  are  used  to  construct  an  estimate  of  multipath- Induced 
Intersymbol  Interference.  The  reconstructed  baseband  waveform  Is  then  delayed 
and  weighted  according  to  the  current  multipath  parameters,  and  this  simulated 
echo  Is  subtracted  from  the  Incoming  baseband  waveform.  The  distinguishing 
feature  of  our  approach  Is  that  an  explicit  model  of  multipath  parameters 
replaces  the  tranvsersal  equalizer  studied  previously. 

An  Important  counterpart  of  the  filter  design  prolem  Is  the  problelm  of 
modeling  non-statlonary  signals.  In  Appendix  E  we  consider  the  problem  of 
estimating  sinusoidal  or  narrowband  signals  with  a  time-varying  center 
frequency.  The  signal  parameters  are  estimated  by  fitting  an  autoregressive 
model  with  time-varying  coefficients  to  the  data.  The  overdetermined  modified 
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Yule-Walker  equations  are  used  to  estimate  a  set  of  constant  model 
parameters.  Some  numerical  examples  Illustrating  the  behavior  of  the 
estimator  are  presented,  and  Its  accuracy  aspects  are  briefly  discussed.  This 
particular  studty  was  performed  only  In  small  part  under  the  current  contract, 
but  was  Included  for  completeness. 

The  work  described  here  represents  an  Important  step  towards 
accomplishing  the  difficult  tasks  of  modeling  and  estimation  of  complex  non¬ 
stationary  signals.  Further  research  Is  needed  Into  the  fundamental 
properties  of  time- varying  processes  and  Into  the  development  of  digital 
signal  processing  techniques  for  handling  such  processes. 


Global  Convergence  of  the  Constant  Modulus  Algorithm 

Julius  0.  Smith 
Bcnjamia  Fnedludcr 

Systems  Control  Technology  lac. 

180t  Page  Mill  Rd.,  Palo  Alto  CA,  94303 


Abstract 

This  paper  proves  global  convergence  of  the  Constant-Modulus 
Algorithm  (CMA)  for  the  case  of  a  real  channel  when  the  model 
order  is  equal  to  or  greater  than  that  of  the  channel  (the  so> 
called  'model-complete*  easel.  The  analysis  is  based  on  an  exact 
fourth-order  Taylor  series  representation  for  the  cost  function 
minimised  asymptotic  ally  by  the  CMA. 

1.  Introduction 

The  CMA  (3.5,6,T.8|  adaptively  equalises  constant- mod¬ 
al  us  communications  signals  such  as  frequency- modulation 
( FM)  and  phase- modulation  (PM)  fortnnta.  The  CMA  adap¬ 
tively  minimizes  a  measure  of  the  amplitude-envelope  dis¬ 
tortion,  inch  as  that  caused  by  multipath  propagation,  us¬ 
ing  some  form  of  gradient  descent  [23.4.9.11|.  The  amplitude 
distortion  measure  is  typically  a  weighted  time-average  of 
an  error  of  the  form  c?  *■  (1  —  where  rit„  it  the 

modulus  of  the  equalized  output  signal  at  time  t.  By  chang¬ 
ing  the  equalizer  parameters  based  on  the  gradient  of  this 
error  measure  at  each  time  sample,  the  channel-induced 
distortions  can  be  eliminated  in  many  situations  (3,8|. 

This  paper  examines  the  convergence  behavior  of  the 
CMA  based  on  gradient  descent.  Three  forma  of  the  CMA 
are  considered,  corresponding  to  three  error  measures.  The 
Snt  is  the  standard  CMA  for  complex  data  [3|,  the  second 
is  a  novel  real-only  form  [9,11|,  and  the  third  is  the  pre¬ 
existing  real-only  form  (fi.8|.  These  will  be  referred  to  as 
cases  1.2,  and  3.  respectively.  Various  extensions  of  the 
CMA  discussed  in  [9,11|  are  incorporated  also. 

Section  2  gives  the  CMA  problem  formulation,  and  sec¬ 
tions  3  and  4  rwpeetively  address  asymptotic  bias  and  con¬ 
vergence  for  gradient-based  implementations  of  the  CMA. 


aal  Model 


Let  yi  denote  the  received  signal  for  f  »  1, 2, . . . ,  T.  In 
the  complex  case,  yt  is  assumed  to  be  of  the  form 

yt  *  fftf(d)* i  +  H9t[d) «i  +■  (l) 

This  work  was  supported  by  the  U.S.  Air  Force  (AFSC),  Rome 
Air  Development  Center,  under  contract  no.  F30H02-84-C-0016. 


where  =»  m, is  the  transmitted  signal  (ot  is  the 
real-valued  information-bearing  signal),  ut  is  additive  white 
noise  at  the  channel  input,  m  is  an  interference  signal 
“template*  (assumed  known),  and  H^d)  and  £f«y(d)  are 
the  unknown  linear  time-invariant  channel  filters  associated 
with  xi  and  a«  respectively: 


•  ■e-f  * 

.4(d)  ^ a(d  +  <ijd*  +  '••  +  (2) 

8|f)  ^  hid  +  bjd*  +  •••  +  4n*d"* 

C[d)  ^  1  +■  etd  +  e2d*  ■“•••+  c«.d^ 

where  are  real.  The  unit-delay  operator  d  is 

defined  by  d*lt  —  r(_a  for  an  arbitrary  signal  tt.  We 
assume  the  modulus  mI(  at  |*,j  is  either  constant  or  known 
for  all  t. 

In  the  'real-only*  caae.  the  transmitted  signal  is  assumed 
to  be  of  the  form  r(  a*  m„  coe(#t).  and  ut  and  u,  are  real 
interference  and  additive  noise,  respectively.  We  assume 
Oe  »  wef  +  Vt  where 

\dvt/dt\<.*t  (3) 

so  that  positive  and  negative  frequency  components  arc  not 
aliased  together. 

Error  Criteria 

The  CMA  minimizes 


,  T 


with  respect  to  the  equalizer  parameter-vector  i.  where 


|  *A8 )  |  -  m;,,  Case  1:  Complex  Signals 

1  i  «-i 

•'{i)  "  2 1  ^  E  *?-.<*>  -  <■  «■»  —  Signals 

if(d)—  Case  3:  Simplified  Reai 
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is  the  instantaneous  error  for  the  three  catee  being  con¬ 
sidered.  and 


*«(f)  A  #^d>[w -£.,(<*)«« 

H  ld\  A  £lii  it  ld\  A  *0£ 


•4(^1  Aii.  -f  4| d  +  iji*  +  •'•  + 

fltd)A  M  +  M*  +  -"  +  6*,<f*‘ 

C(d)  A  1  +  (i^  +  ctd*  +  •  •  •  + 

d  ^j^ao.iii . a*..  it . i*,,fft . 

fo  4[««<«t . «*..*» . **».<! . 
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I<bbi  3  (Blaa).  For  any  gradient-based  CMA.  cases 
1  and  2  are  locally  asymptotically  unbiased  in  the  model- 
complete  ease,  the  minimum  error  being  /(f0)  »  0.  The 
Simplified  Real  CMA  (case  3)  is  generally  asymptotically 
biased. 

Proof.  Let  /(f)  ^  limy— «  £/(£  T),  where  J[S\  T)  a 
given  by  (4).  (The  factor  of  1/2  is  introduced  to  simplify 
later  expressions.)  .An  unbiased  minimizer  can  be  obtained 
only  if  /*(f*)  ■»  0.  We  need  to  show  that  J’(90)  .  0  in 
eases  1  and  2,  and  that  /*(fo)  ^  0  in  case  3.  We  have 

/(f)  A  As)  A  £.*«< fK<f )  (8) 


where  <((f)  is  given  by  ($),  and  e^f)  A  Jiiirj),  where 


If  <*«  >  n€,  then  *«'+*  A  OVfc  >  0.  and  similarly  for 
*»,Ae-  la  (the  real)  cases  2  and  3.  the  modulus  m„  is  re¬ 
placed  by  the  root-mean-square  srTt  (assumed  constant  or 
known).  The  sequence  »«  provides  a  non-negative  weight 
function:  typically.  wt/wt+i  m  \  —  constant  for  a  fixed  ex¬ 
ponential  decay  of  past  data  influence.  The  signals  it,  y«,  u, 
are  complex  in  case  1.  and  real  in  eases  2  and  3.  while  the 
equalizer  parameters  a,-,  and  c,  are  real  in  every  esse. 

Case  2  sums  the  '‘instantaneous  error"  i,  over  the  ■‘ins¬ 
tantaneous  period'  P,.  This  error  definition  for  the  real- 
signal  case  is  essentially  the  same  as  the  complex-signal 
caae  if  the  imaginary  part  is  constructed  at  the  receiver  by 
delaying  the  real  part  oac  quarter  of  the  instantaneous  per¬ 
iod.  Measurement  of  the  instantaneous  period  is  feasible  as 
long  as  (3)  holds.  We  assume  that  the  sum  of  signal  samples 
x,  from  time  t  to  1  —  P,  •+• 1  is  exactly  zero:  in  practice,  oae 
may  wish  to  interpolate  it  to  achieve  this  [13|. 

Case  3  is  the  previously  studied  real-only  algorithm  [4.3|, 
and  here  it  is  viewed  ss  an  approximation  to  case  2.  The 
simplification  introduces  asymptotic  bias  and  a  changed 
asymptotic  parameter  variance.  The  next  section  contains 
formulas  which  can  be  used  to  compute  the  asymptotic  bias 
and  variance  (the  variance  being  inversely  proportional  to 
the  curvature  of  the  error  surface  at  the  minimum  point 
!4)). 

3.  Asymptotic  Bias 

Definition  1.  The  weighted  tim^average  expectation 
operator  is  defined  by 


where  w«  is  a  fixed  ooo-aegative  weight  function. 

Definition  2.  The  nsedcf-eempfctc  case  of  the  CMA 
is  defined  ss  the  ease  n.  >  it.,  ri4  >  n4,  n«  >  ne,  and 
i/t  ■  0.  That  is.  the  model  coefficients  can  be  set  to  exactly 
represent  the  true  channel  and  interference  given  ’jt  and  u». 


Case  1 

ftxiit)  A 

i  p'-1 

p  it-pi/t-p. 

Case  2 

it  it, 

Caae  3 

and 

it  —  £ff.  i 

J  -  J  J  . 

l  *M  •  r-t  — ' 

*/C(d) 

fit  «■'  (fit.  Jr— i  •  •  ■ 

-i<-i . 

Note  that  g  is.  a  linear  operator. 

Now.  /'(f)  *•  0  itf  CtcrifM  it^f  f)  »  0.  ic  cases  l  and 
2.  «.(fo)  a  0  -  A*,,  "*  0. 

In  case  3,  we  have 

/'(fo)  —  E«[if(fc) 

It  suffices  to  show  /'(f#)  yt  0  for  constant  <rx,  =“  <rt, 
no  interference,  and  degenerate  channel  f.  »  1.  In  this 
ease,  ft  a  x,  m  m»co«(o,).  mx  at  \/2<r*.  and  Si  = 
z«.  The  error-surface  gradient  is  then  /'(f«i  =*  E,x}  - 
Setting  m,  —  1  and  Pt  =  we  find  =« 
E«cos3(wef)  =■  1/2  while  Etxf  »  E»cos4(wcD  =»  3/8  A 
f  *  (1/2)*.  Thus,  case  3  necessarily  yields  a  biased 
solution  under  gradient  descent. 

4.  .Asymptotic  Convergence 

Definition  4.  An  m-dimensiooal  matrix  R  of  order  n 
is  defined  as  any  scalar  function  over  a  set  of  m-tupies 
/?(•" i ,  ij, ....  ini  where  each  index  i,  ranges  from  1  to  n. 
j  ■<  1, 2,....m.  Such  a  matrix  will  be  called  an  (m. n| 
matrix. 

Definition  3.  The  /-product  of  an  (m.  n)  matrix  R  times 
a  (1,  n}- matrix  x  (n-vector)  is  defined  ss  the  |m  —  1.  n|- 


8 


matrix 


^  »t . !<><«  I12) 

•7-‘ 

Definition  9.  Aa  (m,  n)  matrix  R  is  said  to  b«  norum- 
gnlar  if  its  /-product  with  any  nonzero  «»vector  is  no o zero 
for  ail  /,  i.e..  RJ*j  —  0  for  soms  j  implies  t «  0. 

Definition  7.  A  signal  y»  is  said  to  b«  penietently 
exciting  (PE)  of  order  (ns,  is)  if  the  (ns,  n  (-matrix 

Rf  A  .  *>  ™  l . »  1 13> 

is  aonsinguiar.  This  definition  can  be  regarded  as  apply¬ 
ing  to  the  received  realization  jr«  of  an  underlying  random 
process.  Normally  we  expect  such  a  random  process  to  be 
PE  with  probability  one  (wpl).  Note  also  that  the  weight¬ 
ing  »(  used  in  the  time-averaging  operator  St  can  affect 
whether  yt  is  PE. 

Definition  9.  A  signal  yt  is  said  to  be  peniotently 
exciting  of  ordtr  (4,  is)  w tlh  respect  to  the  scalar  function  g 
if  the  (4.n(-matrix 


/2J{». *.  11  £  E*g(9t-,gt-j)9(si-kyi-t}  ( M) 

is  nonsingular. 


and  define  g{  • )  as  in  (9).  Expanding  /(0)  in  a  Taylor  series 
about  0*  gives  (exactly) 


m  »  /(0*)  +  0  +  iC,  **(*?*,)  +  et(0*)p(r?) 


jE«3<iJit)jrtx?)  +  ) 


(17) 


and 

A*)  “  S«^*t’ii)j(*?i5e)  +  ei(0*)j(irPt) 

+  j;Ei0<*i  )S<*tVi)  + 

+  (18) 

-  jEs[(*  -  •')*  +  0’)][irtiSi*f  )0 

+ Stp^«t(0*)^tv?rj(0  -  0*) 


Consider  first  cases  1  and  2.  Without  loss  of  generality  we 
can  set  0*  ee  0«  in  (18).  This  gives  immediately  /'(0O)  ” 
/,(—0t)  tm  /'(0)  —  0.  Thus  we  have  found  three  stationary 
points  for  eases  1  and  2.  We  now  show  that  these  are  the 
only  stationary  points.  The  t th  row  of  equation  (18)  with 
0  set  to  0a  (valid  for  cases  1  and  2  only)  can  be  rewritten 


Definition  9.  A  convergent  gradient  descent  algorithm 
is  any  iterative  algorithm  for  0  (of  the  form  0k+t  »■  /*(0fc)) 
[21  which  converges  to  either  a  stationary  point  0*  of  the 
error  surface  /(0)  (in  which  case  J*(0*)  *■  0),  or  to  a  point 
on  the  constraint  boundary  for  0  (if  any).  Normally  this 
property  is  obtained  by  using  a  diminishing  non-summable 
step-size  (such  as  l/k)  in  the  gradient  descent  iteration  [4|. 

Theorem  10  (Global  convergence).  In  the  model- 
complete  case,  yt  persistently  exciting  of  order  (4s«)  with 
respect  to  g  (of  (9))  wpl,  no  interference,  it*  m  ne  »  0. 
then  any  unconstrained  convergent  gradient  descent  CMA 
will  converge  with  probability  one  to 

0*  —  ±*»  »  ±(*».  at ,  a*, ... ,  a*. ,  0, . . . ,  0|r  (15) 

in  cases  1  and  2,  for  any  nonzero  initial  parameter  vector 
0(0)  yd  0. 

Proof.  The  stationary  points  of  a  gradient  descent  al¬ 
gorithm  occur  at  points  0  where  the  gradient  J'{9)  is  zero. 
The  first  goal  is  to  find  ail  such  points. 

Let  0*  denote  any  local  minimizer  of  the  cost  function 
J{t)  given  by  equation  (8).  By  definition,  7’(0")  as  0.  By 
lemma  3,  0  *  0«  is  one  such  minimizer.  Let 

h  Asfrf  =  do) 


Aim  -  ;  E  E  £  i.  M  ( 

*  I— »z— I  Jr— l 

where  a,-  »  (0  +  0e)[«1,  3j  —  0(/|,  and  tp  —  (0  -  0O)W- 
The  (4,  ri. (-matrix  flJi'i./.M  ^  E«j< )gi Vt-k’Jt-t ) 
can  be  interpreted  as  a  type  of  four-dimensional  covariance 
matrix.  We  see  that  a  sufficient  condition  for  the  three 
solutions  0  €  (9<  0o.  — 0o)  to  be  the  only  solutions  is  to  have 
R*  be  nonsingular.  But  this  holds  whenever  yt  is  persis¬ 
tently  exciting  of  order  (4.n«)  with  respect  to  g.  Hence. 
O,0g,  and  —0e  are  the  only  stationary  points  in  cases  1  and 
2  with  probability  one. 

It  remains  to  be  shown  that  0  —  0  is  an  unstable  station¬ 
ary  point,  while  0  —  ±0*  are  stable  stationary  points  tor 
the  gradient.  A  stationary  point  0  is  stable  if  the  matrix 
/*(0")  is  positive  definite,  and  unstable  if  /"(0  )  is  negative 
definite  [2|.  It  is  straightforward  to  show  that 

A(0)  —  -y"(±0o)/2  (20) 

That  is,  the  curvature  of  the  error  surface  at  9  »  0  is  equal 
to  —1/2  the  curvature  at  0  «  Jo,  and  the  curvature  at 
—0«  equals  that  at  0#.  It  is  always  true  that  /"(0O)  >  0 
(see  the  third  term  in  ( 17)),  and  yt  persistently  exciting  of 
order  (4,  n«)  with  respect  to  g  implies  9q )  >  0  (positive 
definite). 


Thus  /*|0)  <  0.  making  9  »  o  an  unstable  stationary 
point,  jin  ail  current  forms  of  the  CMA.  9  can  never  be 
allowed  to  equal  zero,  for  this  will  in  fact  “freeze'  a  gradient 
descent  algorithm.) 

The  two  remaining  solutions  9  <m  ±09  are  $tatU  station* 
ary  points.  (/*(— do)  ■»  J"{9 !»)  >  Oj.  Under  the  assumed 
conditions,  9  ”  ±9q  are  the  onlf  stable  stationary  points 
of  the  complex  or  the  real  CMA  (cases  1  and  2).  Since 
in  practice  bf„(d)  is  divided  through  by  d«,  the  sign  am¬ 
biguity  in  9"  becomes  a  sign  ambiguity  in  H t*U);  a  sign 
change  in  a  signal  is  normally  negligible.  It  is  interesting 
to  note  that  this  sign  ambiguity  corresponds  to  the  phase 
ambiguity  in  the  complex-channel  case  (3|. 

Further  Enecsions 

Theorem  10  immediately  extends  to  the  case  which  in¬ 
cludes  interference  canceling  (a*  >  n*  >  0).  Similar  global 
convergence  results  may  hold  for  the  model-complete  cases 
in  which  nc  >  0;  for  example,  the  techniques  in  [1|  might 
go  through.  It  has  been  shown  (fi|  that  either  noise  (iq 
0)  or  model-incomplete  interference  (Hap  yd  //,,  VP)  will 
cause  bias  in  the  parameters.  For  case  3.  theorem  10  can  be 
extended  to  prove  global  convergence  to  a  biased  solution 
in  which  the  bias  can  be  simply  approximated  u  in  [6|. 

Convergence  results  for  specific  algorithms  can  be  ob¬ 
tained  using  the  analytical  approach  described  in  [1|.  In 
the  model  incomplete  case,  the  general  model  (3)  should 

converge  to  a  local  minimum  of  the  error  surface  J(9):  the 
number  of  sub-optimal  local  minima  can  be  large  in  the 
model-incomplete  case. 

We  expect  that  if  3,  am  and  <1  is  randomly 

distributed  with  almost  any  non-discrete  distribution,  then 
Ry  will  be  nonsingular  with  probability  one  for  any  number 
n  of  parameters.  It  seems,  however,  that  for  good  numeri¬ 
cal  conditioning,  further  restrictions  are  necessary  on  the 
modulating  signal.  For  example,  it  might  be  appropriate 
to  require  jt  to  possesa  at  least  a  distinct  frequencies  of 
high  spectral  power,  analogous  to  the  situation  in  least- 
sqnara  system  identification  (4|.  Further  work  is  necessary 
to  specify  precisely  the  modulation  characteristics  which 
maximize  the  equalization  accuracy. 

The  error  minimized  in  the  model-compiete  case  is  ex¬ 
actly  described  by  a  fourth-order  Taylor  series  (cf.  equation 
(17)).  A  fourth-order  type  of  gradient-descent  algorithm, 
analogous  to  Newton  s  method  for  least  squares  problems, 
should  yield  the  fastest-converging  algorithms.  To  this  end. 
note  that  all  solutions  to  the  ensuing  third-order  "vector 
polynomial’'  for  the  gradient  can  be  expressed  in  cloeed 
form. 

ConclMtow 

Some  convergence  properties  of  the  Constant  Modulus 
Algorithm  (CMA)  for  channel  equalization  were  described. 
Subject  to  mild  restrictions  on  the  modulation  signal,  the 


unbiased  versions  of  the  CMA  can  only  converge  to  plus  or 
minus  the  true  solution  in  the  model-complete  case. 
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Abrtwrt 

The  Coaataat-Moduiua  Algorithm  (CMA)  computes  ted 
applies  sa  adaptive  channel  equalizer  for  constant- amplitude 
death  such  ss  frequency-  sad  phas^  modulation.  This 
paper  prostata  several  ccteuimas  to  the  CMA  including  OR 
equalisation.  s  real-signal  ▼tnwa  having  praportiss  as  good 
as  tht  comp  its  version.  use  of  the  Gaase-N'ewtoo  method  ia 
plan  of  gradient  donsat.  latarfaroan  rejection.  sad  mart. 

Ths  Coascaat  Mod  alas  Algorithm  (CMA)  was  introduced 
hr  J.  &  Thdehlor  k.  si.  (4.7.3,9.13)  as  a  mothod  for  adap¬ 
tive  M  -—«■»«"  >yp—  xt "y«i« 

The  CMA  calibrates  a  Uaoar  ehaaaoi  equalizer  hf  seeking 
to  make  ths  equalized  signal  hare  coastaat  mod  alas.  Sash  a 
technique  «n  ho  oaod  with  freqeeacy-modeiatioa  (FM)  tad 
phaewmoduiatma  (PM1  comanaicmiioai  systems.  when  ths 
amplitude  envelope  is  coastaat  or  known  ia  tho  ahaoaeo  of 
diatortioa.  Tho  CMA  ia  oat  example  of  adaptive  ehaaaoi 
equalisation  haard  oa  restoring  invariant  praportiss  of  the 
uadis  tort  od  sigasi  jo). 

Ia  [e|  tho  roc arod  signal  is  assaaud  available  as  a  com¬ 
plex  signal  of  tho  form  s'**.  Ia  this  eaao,  tho  modalaa  is 
exactly  coastaat  whoa  ehaaaoi  distortion  an  i  be  eat.  Ia 
[1.12),  a  real-only  version  of  tho  CMA  ia  devdopod.  Moat 
raeaatly,  tho  CMA  has  been  ext  ended  to  known-amplitude 
(aa  opposed  to  constant- ampiilade),  sad  to  maltiehaaasl 
equalization.  Ia  all  eases,  the  equalizer  is  formed  by  apply- 
iag  a  Sxud-step  gradieat-daMtat  algorithm  similar  to  the 
Wldraw  IMS  algorithm  [ij. 

Tha  purpose  of  this  paper  is  to  praseat  several  extensions 
to  tha  basic  CMA 

e  The  ehaaaei  equalizer  is  aa  arbitrary  rational  liter 
(poles  aad  teres),  or  infinite-impulse- response  (HR) 
liter,  rather  thaa  beiag  eeastraiasd  to  a  3nit»impuls* 
reepoaae  (FIR)  ail-iaro  Altar  u  ia  eariiar  works.  Use 
of  aa  til- pole  equalizer,  for  example,  allows  aa  exact 
inverse  to  FIR  ehaaaei  distomoas  such  ta  multipath. 

0  A  saw  typo  of  real-only  algorithm  is  derived.  The 
sew  real-only  version  is  based  oa  a  more  desirable 

.’.sis  «er*  -as  suooartM  3y  cite  J.S.  Air  Pome  (APsC). 
under  Contract  to.  .*30602-dt-C-M1S. 


error  definition  whieh  is  equivalent  to  tho  complex 
eaao  when  the  iastaataaeons  carrier  period  is  used. 

s  The  Recursive  Gauss-Newton  method  (RGN)  [3.81 
replaces  the  gradient  descent  method  for  obtaining 
tha  equalizer  coefficients.  This  modification  can  im¬ 
prove  the  convergence  rate  and  ssyraptouc  accuracy 
of  the  equalization. 

e  A  new  typo  of  RGN  algorithm  is  presented  [S|  which 
aliowe  a  trade-off  between  properties  of  the  more 
recent  recursive  forms  aad  the  origiaal  'offline'  or 
•batch’  forms. 

e  A  technique  for  removing  partially- known  interference 
ia  incorporated  into  the  algorithm. 

Ia  sactioa  2.  the  basic  problem  formulatma  ia  presented, 
and  sactioa  3  derives  the  algorithm. 


Signal  Model 

Let  pi  denote  the  received  signal  for  I  •  1. 2. ....  T.  la 
the  complex  case,  j,  ia  assumed  to  be  of  the  form 

n  «  tfe*l d)*f  -  H^d\u,  -  H,9{d)v‘  1 1) 

where  si  >■  is  the  traasmitted  signal  (oi  ia  the  real- 

vaiued  information- bearing  signal),  is  sddiuve  noise.  a, 
a  aa  interference  signal  (aaaumed  at  least  partially  known), 
aad  Jftfii)  sad  ffay(d)  are  the  linear  time-invariant  chan¬ 
nel  liters  associated  with  r<  aad  at  respectively: 


«..;«>  i  ?£ 


•■Md)  4«e  *«id -u-d*  - - (3) 

m  £  M-M*  - 

CTd)  ^lvctiv «jd*  -  •  •  •  — 

where  *,)  are  real.  The  unit-delay  operator  d  ia 

defined  by  the  rciation  d"n  ■»  *«_*  for  aa  arbitrary  signal 
>i.  We  assame  the  modulus  oiv,  of  the  transmitted  signal 
*i  is  known  for  ail  1. 

The  novel  features  introduced  so  far  ar  the  provision  for 
tfapid)  aad  ««  to  model  partiailr-cnewu  interference,  aad 
owe  Air  Cevelocnant  Center . 


rnevious  root 

is  BLANK 
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the  polymaul  C\  d)  ia  the  channel  transfer  faaetioa  H^d) 
tor  mod  disc  FIR  distortion  typos. 

Ia  tbs  east  of  mi  signals.  tbs  transmitted  signal  is  is- 
sumsd  to  bo  of  tbs  form  <i  «■  mt,  eosf  #t)  sad  «»  sad  vi  in 
rsal  iatsrfsrsaes  sod  additive  notes.  respectively.  Ws  also 
sssaas  Si  whwe  \9*i /dt\  <  w*.  Ia  this  csss 

tbs  known- modulus  proporty  ess  bs  rtplsesd  by  tbs  knows- 
eaeefsp*  propsrty.  Homwr.  to  rsuia  a  lent  spurts  prob¬ 
lem.  ws  oss  isstssd  aa  cypress  rasfcfy  buss  cssrsp *  power 
propsrty.  Thus,  ia  tbs  rssi-sigasl  csss. 


(3) 


wbsrs  /*<  is  tbs  ‘iastaatsnsous  period’  of  *<,  aad  <r^t  is 
dsdasd  as  tbs  ‘iastaotaasous  powsr,*  or  'intift  power* 
of  S|.  For  simpikity.  ws  Latsr  ass  equality  ia  (3).  Tbs 
perturbaooas  ia  (3)  dus  to  disersts  Pi  aad  moduiatioa  by 
«t  esa  bs  compensated  if  dssirsd  by  waveform  iatarpolatioa 
aad  rsdefinitioa  of  cootunous  P,  to  absorb  tbs  earner-cycle 
rapport  distortioa  oadsr  integration. 

Tb»  £aat 

Bolow  ws  dadas  tbrts  am  entsru  for  tbs  CMA.  Tbs 
tbrss  mar  ehtsris  correspond  to  tbs  csss  of  complex  s, 
(4).  rsal  2i.  aad  aaunpiifiod  ▼stsm*  for  rsii  z(  [8,121.  These 
will  bs  rsfsrrsd  to  a  rasss  1 sad  3,  respectively. 

Cass  l:  Error  for  Camolss  Signals 
la  tbs  eoapiss  csss.  tbs  CMA  minimum 

r 


A  ri  A  5".  •*««?.  j(*j 


(-*» 


with  impset  to  tbs  equaiixor  parameter-vector  #,  wbsrs 


i«.»(«)4p.(f)|  -ml, 


<?WJ  r,  s 


ft  -Id\  i  Zll! 

****** 


(S) 


.4fd)  |^«s  bid*  ■+••••  * 

3(dl  ^  M  M*  *> •••  ■* is,d"> 

£(«fl  4  l+e,d* ijd3  l^d"‘ 

*  Aj^ds.  A . ds„»t . S*A.---.*«.| 

aad  «i  is  a  aon-aegative  weighting  function.  Tbs  signals 
st.  n, ««  art  eoapiss,  wbils  tbs  equaliser  parsmstsrs  i,, 
aad  e*  ars  rsal. 

Ws  tbow  ia  (u;  that  if  m  0.  n»  »  x«  m  9,  nt  > 
m«.  «j  >  nv,  sad  ii  is  ptniiUntlj  eeiiitty  ia  s  hovel 


sense.  tbsa  tbs  oaiy  psrsmstsrs  9  wbieb  locally  minimise 
/,(#;  n  ars  4  m  ±ds  (wbsrs  4s  is  extended  with  terra  to 
tbs  siss  of  4  if  necessary).  This  is  valuable  to  know  trace 
it  asnaa  that  tbs  only  minimises  of  A(4;  T)  ars  global 
minimum.  Ia  aneb  a  esss.  gradient  descent  methods  esa 
bs  used  with  confidence  ia  tbs  so-called  -model- complete’ 
esss  (La..  when  tbs  ebanosi  model  esa  exactly  describe 
tbs  true  ebaaasl  distortioa).  Furthermore,  the  psrsistsnce- 
of-cxej tattoo  condition  bolds  for  ‘almost  all"  information- 
I  idea  modulation  signals  ot. 

Cass  2:  Error  for  Rsal  Signals 

Ws  dcdai  tbs  rssi-ooiy  case  as  minimizing 


wbsrs 


.  A- • 

**.t iilln-  £  i?-,(4)-oi, 
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aad  all  other  quantities  are  formally  tbs  same  as  ia  (a)  with 
A,  yi,  *<  real.  Thus  tbs  maw  difference  between  the  com¬ 
plex  aad  real-only  cases  is  that  tbs  ‘‘instantaneous  error* 
i|.i  is  replaced  by  aa  error  ej.t  wbieb  is  aa  average  over 
tbs  ‘instantaneous  period’  Pi.  This  definition  for  tbs  real- 
li goal  csss  is  cqsrcataif  ia  principle  to  tbs  complex-signal 
ease  when  tbs  imaginary  part  is  constructed  it  tbs  receiver 
from  tbs  rsal  part. 

While  (71  may  mm  like  tbs  natural  error  to  use  in  the 
real-only  ease,  it  requires  aa  wtimnte  of  tbs  instantaneous 
period  Pi  wbieb  depends  on  tbs  modulation  signal  jt.  Tbs 
counterpart  ia  tbs  complex  CMA  is  to  define  tbs  imagi¬ 
nary  pan  as  a  quarter-cycle  deiay  of  .be  real  part,  where 
tbs  quarter-cycle  delay  is  ons-fourtb  of  tbs  current  rastaa- 
taasoua  period.  Note  that  measurement  of  the  instan¬ 
taneous  psriod  requires  using  a  earner  frequency  aueb  higher 
tbsa  tbs  signal  bandwidth. 

Aa  altera atrvs  to  tracking  tbs  instantaneous  psriod  is  to 
Sx  Pi  to  a  value  wbieb  spaas  many  eyeies  of  tbs  carrier 
under  ail  moduiatioa  conditions.  Tbs  counterpart  in  tbs 
complex  CMA  is  to  use  tbs  Hilbert  transform  to  create 
tbs  imaginary  pvt  of  r«  ( as  ia  typically  doas  in  practice*. 
Averaging  over  many  earner  cycles  has  tbs  advantage  of 
aot  requiring  aa  estimate  of  instantaneous  period,  hut  the 
disadvantage  of  limiting  tbs  potential  convergence  rats  of 
tbs  eqoaiixv  parameters  4. 

Casa  t  Simplified  Error  for  Rsal  Signals 

A  simplified  real-only  CMA  [3.12|  minimises 


.  T 

Ai*  HaeT  9) 

*  e— i 
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where 


iu(#i  io> 

While  this  ia  aot  obviously  u  appropriate  error  criterion 
for  non-consunt  signal*  **,  it  cm  b«  shown  (U|  that  (9)  ha* 
asymptotic  properties  r*r y  similar  to  thcae  of  (7).  Intui¬ 
tively,  the  roasM  is  that  tbs  periodic  fluctuations  of  the 
error  is  (9)  at  the  ioeuattaeose  earner  frequency  average 
oat.  The  error  limplifleatio*  uitrodace*  asymptotic  bias 
(oftem  appraximable  by  a  topic  reals  factor  [8.X1|]  aad  a 
riifbtly  different  asymptotic  parsmster  variance  relative  to 
tbs  optimal  versions  (41  aad  (6)  (ll|. 

3.  The  .■Vi  corn  h  to 

Prior  trescmeacs  of  tbs  CMA  bare  bees  baasd  oa  pradi- 
tat  descant  [4.7.9,9.13].  G radical  deeeeat  cm  be  acea  11 
a  linearisation  metbod,  while  Newtoo  d  cm  eat  mores  each 
itsratioa  to  tbe  bottom  of  a  quadratic  local  approximation. 
Ia  Newtoa  descant.  close  to  tbe  aiaiomm  point,  coover- 
peace  is  quadratic  while  (radical  dmeeat  converges  linearly 
[3|.  Taos.  Newtons  metbod  converges  faster  aesr  a  local 
miaimam  (aseomiap  tbe  ITmsiai  is  aousero). 

Tbe  Grsts-Vewtoa  metbod  (3)  is  a  robast  approximate 
Newton's  metbod.  It  is  preferable  orer  Newtoa  *  method 
when  time-recursive  algorithms  are  desired.  Ia  this  sec. 
tioa.  we  derive  the  Gauaa-Newtoe  metbod  as  applied  to  tbs 
CMA. 

Tbs  CMA  coat  fuactioa  cm  be  espreaaed  a 

■  .  r 

A*  rj  A  57  £*»#•)  (i°> 

I*  I 

where  w«  is  a  non-negative  real-raised  weighting  function, 
sad 

|  dc(f)  j*  -  mj,.  Cass  U  Complex  Signals 

hi*)  -  5  i  £  -#*„  Case  3:  Real  Sipaals 

‘  »-* 

*t(  j)  —  #i,,  Cass  X  Simpliffed  Real 

(U) 

Tbs  facton  of  1/3  ia  ( 10)  aad  ( 11)  bare  beaa  introducad  to 
simplify  later  erpraaaioaa.  Let 


4<«)  4 


(12) 


Tbe  Casae-Newtoa  1GN)  Metbod 

Let  f*  desote  a  local  miaimiaar  of  /( £  D.  Whan  the 
trror  «<  cm  be  dricts  to  taro  as  '$  approaches  #*.  or 


aad  <1  ( 9*1  are  uncorreiated.  tbs  Hessian  cm  be  dosety  *p- 
proximal ed  by  Tba  Co ssa-Newisn  [3|  method  ia 

detaad  accordingly  as  follows.  Gi«ea  u  iaitial  parameter 
estimate  ?«,  carry  out  tbe  foilowiap  iteratioa  until  coarer- 
peace  ia  achieved: 

.  T 

Rf9.i  T)  -  j  £  "JMiW.f  *  '*(«»  T)  (X3J 
9,v.  —  «.  -  «“*(#.;  T)A*i.T) 


The  IT  win  a  approximation  serres  two  purposes.  First,  it 
the  seed  for  sacond-ordar  differentiation— tie 
pradient  alone  drive*  tbe  parameter  updates.  Second,  tbe 
matrix  R I9;;D  ia  more  likely  to  be  positive  del  nice  and 
inreruble  thM  is  J"[ih.  T)  [31. 

Exact  Raenrsirt  GN  Updates 

Tba  Gauaa-Newtoo  algorithm  (13)  cm  made  recursire. 
where  the  parameter  animat*  i  is  updated  for  each  t.  We 
now  consider  only  one  peas  through  the  T  data  points. 
Therefore,  tba  pasa- number  •  in  (13)  ia  dropped  for  aota- 
lioasl  simplicity  Tbe  initial  value  ♦,  from  tbe  previous 
pass  ia  denoted  f«.  and  the  Inal  minute  d,„,  obtained  it 
the  end  of  pas*  ■  sow  becomes  #7- 
Oedoe 

r  ,  r 

Hr  *i«t(9eli't(9e)r  gT  A  —  Y'  wr«'ii9o)c«i9o i 

*T  *T 

I  HI 

Then 

R.-X,A.t*«'t(9,)i'(9,)r  , 

Gi  ■■  XtGi.|  *■ 

when  X|  A  wi_i/wi  is  known  as  the  forgetting  factor." 
Typically  X,  is  lead  between  0  and  1  to  obtain  a  flxed  ex¬ 
ponential  fading  of  poet  observation*  (with  an  approximate 
time  constMt  of  1/(1  —  X)  tampion.  It  cm  be  shows  ;5j 
that  for  A  any  integer  between  1  and  t. 


!i.l9el**,.(9eir(fr-a-9e)! 
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Tbia  recursive  form  is  exactly  equivalent  to  :be  'off-line* 
counterpart  (13). 

Block- Recursive  Gauae-Newton  i  BRGN) 

Tba  BRGN  [5|  ia  obtained  by  raplacinp  tba  initial  con¬ 
dition  #0  by  recently  computed  estimate*  j,  u  periodic  in¬ 
tervals.  Whan  ie  is  repiacad  by  we  effectrveiy  split  :be 
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In  ail  earn 


data  into  more  than  om  "batch.*  The  lot  batch  consists 
of  data  tip  through  time  t,  and  it  ia  used  to  produce  the 
estimate  The  second  batch  begins  at  time  t  —  1.  and  it 
ia  initialized  to  the  — produced  by  the  dm  batch. 

Let  P,  denote  the  "batch  site*  ia  aa  otf-Uns  Cause-Newton 
method  implemented  recursively  eta  equation  ( IS).  Assume 
for  simplicity  that  A  divides  Pt-  Then  ( 18)  becomes 


it  ■  it—p,  —  RT‘Ct 


[e,i#i.ft)  -p  if  -  #i-ft) 


This  adaptive  algorithm  allows  interpolation  between  the 
properties  of  an  "off-line’  or  "hatch*  optimization  aad  the 
truly  umwrecuniv*  algorithm*  which  are  used  extensively 
in  the  idenuffcation  literature.  The  practical  importance 
of  the  8RGN  for  the  CMA  is  that  it  allows  paraoMter  es¬ 
timates  per-pened  ia  the  rcaf-owfy  case.  The  earner  period 
eaa  be  estimated  by  BRGN  itself  or  it  can  d  served  from  the 
demodulated  signal  ia  a  "daemon  directed"  mode. 

Recursive  Gauss-Newton  IRON) 

3y  setting  the  block  stze  Pi  to  1  (forcing  h  to  1),  the 
BRGN'  reduces  to  the  rersrseee  Geeee-iVewSen  melted  (RCN) 
which  is  somewhat  of  a  standard  method  for  ARMAX  sys¬ 
tem  idenohcation  [8): 


*«  ^  Zt  * 

fit  4 [»••••■  •  »»-«••  . ”*«-«»•  “**-» . 

•  . *"••*« . . *"* 

1*1  we  <fi{  k  fit  -  C,  ,j(_,  -  Crfit-t  - - *"•  r'l-"* 

120) 

These  quantities  determine  the  error  *i  aad  its  gradient  r', 
needed  to  drive  the  Gaus»tNewtoa  algorithm. 

Because  the  recursions  for  it  aad  js{  involve  the  applies- 
tioe  of  the  um^*arymg  filter  l/Cdd).  it  a  accessary  that 
C,id)  be  stable  at  all  tunes.  Our  strategy  for  stahiuty 
projection  is  as  follows:  If  any  root  of  C,U!  lies  on  or  ouuide 
the  drde.  the  roots  are  contracted  uniformly  by  tbe 
factor  |  m  0.04  repeatedly  until  all  roots  are  maids.  Thu 
aagi*«avariaat  projectioo  preserves  the  tuning  of  spectral 
resonances  in  the  equaliser  |e.g.  multipath  auilsl. 

la  the  RGN  algorithm.  4  —  fi,/C[d\,  where  C  a  con¬ 
stantly  being  replaced  by  the  latest  avsilable  estimate.  Since 
i,  ■  computed  using  #t-t.  ths  first  occurrence  of  s  in 
jf/  depends  also  on  #t-i  U  has  been  observed  in  the 
analogous  system  identification  context  that  convergence  a 
accelerated  by  recomputing  i,  using  ft  before  using  it  m 
the  reeursson  for  i{. 

.Algorithm  Summary 


it  -  ie- .  -  )*•<*•- .) 

A  m  X,*,.,  -  K'|i,.,)ef{i,_,)r 

Ljung  [8|  aad  others  have  shown  that,  under  general  condi¬ 
tions.  the  replacement  of  I*  by  even  the  most  recently  trad¬ 
able  parameter  mtimntm  <i.t  does  sot  altar  ths  asymptotic 
convergence  properties  of  the  Gauss-Newton  method. 

To  tailor  GN,  BRGN,  or  RGN  to  a  specific  application, 
the  instantaneous  error  it  aad  its  gradient  if,  with  respect 
to  the  most  recently  computed  parameters  i  art  seeded. 

The  CMA  Error  and  Grad-rot 

The  "Instantaneous  gradient"  of  it  with  respect  to  i  is 
given  by  «^i)  m  i\s,jJx)  where 
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Case  l:  Complex  Signals 
Case  2:  Real  Signals 
Case  &  Simplified  Real 


t{  *  fl  -  ir-ii  Mt£-i - - C|-1  ["iiyf. 
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h«(d)i  -mf,,  Cass  l:  Complex  Signals 
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jjr  £  i1-,li)-<rl,.  Case  i  Real  Signals 
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if(#l  —  e;,,  Case  X  Simplified  Real 
£*{£(£,},  Case  l:  Complex  Signals 

JTt  S  ci-pit-e.  Case  2:  Real  Signals 

iii!,.  Case  3:  Simplified  Real 
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The  Cmmui  Modal oi  Algorithm  (CMA)  for  channel 
equalization  baa  boon  extaadad  ia  seeurai  ways.  The  extea- 
mom  warn  aiaad  primarily  at  improving  convergence  rata 
tad  ramtaace  to  noiso.  The  i  tract  oral  cxtaaaio«  iaeloda 
uae  of  the  rac  antra  Gauao-Newtoa  (RCN)  mat  bod  ia  plan 
of  fradiaat  descent.  a  form  of  RCN  which  interpolates  bo- 
twaaa  hatch  tad  resumes  forma,  aa  iatarfaraaca  caaeeiing 
factlicy.  aad  extsasioa  to  meviag-avuraft  chaaaab. 

Wo  steaded  the  CMA  from  a  gradient  dotcoat  to  a 
Nawtoa  daacaat.  Thu  steads  a  Srst-order  Taiyoc-eerim 
ippramutioa  of  the  error  nrfaca  to  a  seeoad-erder  oaa. 
However.  the  error  miaimisad  ia  thmcaaaia  exactly  described 
by  a  fourth-order  Talyor  Sanaa.  A  fourth-order  typo  of  d* 
mat  algorithm  shooid  yield  faatar  Mtmpatt  Nou  that 
ail  soianooa  to  the  third  ordar  *vuetor  poiyaomial*  for  the 
tradiaot  caa  bo  tcprmaad  ia  eioaad  form. 

Ia  caao  3,  thara  ia  the  paaaibiUty  of  adaptively  tracking 
the  iaataataooooa  period  of  the  carrier  Pt.  Period  tracking 
ia  accompliahad  by  adding  Pt  to  the  parameter  rector  # 
aad  redwtviag  the  error  pdiaot.  If  iho  carrier  freqaeacy 
ia  large,  the  algorithm  caa  emilr  track  Pi-  Aa  a  farther 
refinement,  the  iastaataacooa  period  caa  bo  extended  to 
noa-iategw  ralom  mug  the  techniques  diacoaaad  ia  (10). 

Ia  the  caao  of  maltipath  distortion.  the  ehaaad  model  ia 
of  the  form  l+pdr,  where  typically  r  a  taaay  samp  lea.  U  ia 
jaaacaaaary  to  compute  ail  the  parameters  of  a  Urge-order 
inverse  liter  «haa  traafciag  such  a  chaaaai  model.  The 
taro  ceefitatats  caa  ha  sHminacod  from  the  modal.  A  high- 
quality  method  for  idaptiraiy  traekiag  the  (iatarpoiatad) 
delay  r  ia  this  context  is  daacribad  ia  (10). 

la  digital  uoduiatioa  formats,  the  rwitchiag  traasieata 
may  ha  low.  energy  raiatira  to  the  signal  energy  between 
rwitchiag  instants.  Ia  ordar  to  avoid  iU-eoaditioaiag  due 
to  this  uahaiaaca.  equation  (17)  caa  bo  used  to  skip  over 
iatrvbaud  data.  Ia  other  words,  data  acquimioa  caa  be 
imaged  to  occur  only  ia  a  neighborhood  of  tho  switching 
cam. 

Duo  to  the  tima-ehift  structure  ia  tbo  RCN  algorithms 
we  hare  derived,  there  mist  so-called  ‘fast*  versions  which 
require  order  a  computations  instead  of  nJ  as  ia  the  present 
vensoas.  A  point  of  departure  from  which  such  aa  txteo- 
swu  ia  straightforward  is  grvee  ia  (2). 
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Abstract 


The  set  of  finite-order,  linear,  time-invariant  filters  is  derived  by  sampling  lossless 
propagation  through  a  variable-impedance  medium.  This  leads  to  a  flexible  class  of  Cime- 
vai ryiag  filter  structures,  termed  “Waveguide  Filters"  (WGF)  in  which  signal  power  is 
decoupled  from  changes  in  the  filter  parameters.  These  structures  are  “balanced"  in 
the  sense  that  the  decoupling  between  signal  power  and  time-varying  filter  coefficients 
is  maintained  for  each  individual  section  in  the  structure.  In  addition,  limit  cycles  and 
overflow  oscillations  are  suppressed,  even  in  the  time-varying  case,  when  implemented  with 
"passive"  arithmetic.  We  describe  also  a  method  for  enforcing  exact  losslessness  in  the 
realisation  of  an  arbitrary  digital  filter  in  spite  of  the  presence  of  round-off  errors.  Finally, 
the  WGF  structures  can  be  interconnected  in  series  or  in  parallel  in  a  way  which  does 
not  disturb  the  signal/coefficient  decoupling  or  the  power  balance.  Thus,  the  waveguide 
filters  are  very  useful  for  modeling  physical  systems,  and  the  exactness  of  their  physical 
interpretation  enhances  their  suitability  for  the  time-varying  case.  All  results  are  obtained 
for  the  multi-input/ multi-output  case. 


uw 


1.  Introduction 

Digital  filtering  techniques  have  often  been  derived  from  classical  or  “analog" 
techniques  (29j.  Classical  filter  design  has  its  roots  in  "network  theory”  for  describ¬ 
ing  linear  time-invariant  systems  accessed  by  means  of  "ports”  (11|.  Network  theory 
itself  is  a  body  of  mathematics  built  upon  certain  assumptions  [7,13]  which  become 
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true  in  the  limit  et  low-frequencies  according  to  Maxwell's  equations  for  electromag¬ 
netic  propagation  (0).  Thus,  the  theory  of  filters  grew  originally  out  of  the  scalar 
theory  of  wave  propagation. 

Since  the  emergence  of  digital  techniques,  little  attention  has  been  paid  to 
the  close  correspondence  between  filter  computations  and  physical  law.  In  signal 
processing  applications,  we  normally  approximate  directly  some  desired  transfor¬ 
mation  of  the  signal  spectrum,  and  a  true  physical  modeling  is  irrelevant. 

The  mainstream  of  filtering  applications  has  involved  time-invariant  filters 
which  approximate  an  ideal  amplitude  response  such  as  low-pass,  high-pass,  band¬ 
pass,  or  band-reject  characteristics,  or  which  provide  a  desired  phase  response  such 
as  in  equalizers  for  communications  channels  [29|.  In  the  time-invariant  case,  the 
amplitude  response  and  phase  response  completely  determine  a  linear  filter  [20] . 

For  time-varying  filters,  there  is  no  longer  a  simple  description  in  terms  of 
amplitude  and  phase  response.  (A  frequency  response  requires  time-invariance.)  In 
many  cases,  time-varying  filters  have  been  developed  in  an  ad  hoc  manner,  being 
regarded  as  “quasi-static"  in  most  cases.  Such  extensions  require  the  assumption 
that  the  filter  coefficients  vary  slowly  relative  to  the  impulse- response  duration  of 
the  filter.  When  the  coefficients  change  too  rapidly,  unnatural  artifacts  can  occur 
due  to  the  incompatibility  between  the  filter  state  (a  function  of  all  prior  time  in  a 
recursive  filter)  and  the  new  filter  coefficients. 

This  Paper 

This  paper  presents  a  class  of  recursive  digital  filters  designed  specifically  to 
have  the  best  possible  behavior  under  time-varying  conditions  in  the  presence  of 
round-off  error.  We  call  them  “Waveguide  Filters"  (WGF),  because  they  can  be 
interpreted  as  networks  of  intersecting  waveguides.  The  WGF  structures  are  closely 
related  to  the  “wave  digital  filters"  developed  principally  by  Fettweiss  [17,28),  the 
lattice  filter  structures  arising  in  geoscience  and  speech  modeling  [45,31],  and  the 
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“normalized  ladder  filter"  discussed  by  Gray  (30,391.  Waveguide  filters  have  the 
following  characteristics: 

•  The  correspondence  to  physical  wave-propagation  systems  is  exact  even  though 
time  is  discrete.  No  bilinear  transformation  is  necessary  to  connect  digital 
quantities  with  physical  quantities  as  is  usual  in  the  wave  digital  filter  (WDF) 
context  [17].  This  allows  a  priori  choice  of  filter  structure  to  obtain  precise 
models  for  physical  processes. 

•  The  instantaneous  power  at  each  internal  filter  section  is  invariant  with 
respect  to  filter  coefficient  variation. 

•  Generalized  versions  of  the  “Normalized  Ladder,"  “One-Multiplier  Lattice," 
and  other  ladder/lattice  filters  are  derived,  all  having  invariant  instantaneous 
power  in  the  time-varying  case. 

•  The  structures  can  be  coupled  at  a  junction,  cascaded,  loop«*d,  or  branched, 
to  any  degree  of  network  complexity,  and  the  desirable  properties  such  as 
stability  and  power  decoupling  are  retained. 

•  A  synthesis  procedure  exists  for  computing  all-pole  or  pole-zero  sections. 

•  There  is  an  identification  method  for  determining  the  coefficients  of  the 
structure  from  measured  input/output  data.  Similarly,  there  is  a  “linear 
prediction"  modeling  technique  for  these  structures  which  provides  ARMA 
models  for  time  series. 

•  No  overflow  oscillations  can  occur,  even  in  the  time-varying  case. 

a  No  limit  cycles  (also  called  “granularity  oscillations")  can  occur  if  one  of 
many  “passive"  numerical  round-off  strategies  is  employed,  even  in  the  time- 
varying  case.  In  the  simplest  case,  the  passive  round-off  strategy  reduces  to 
magnitude  truncation  (or  truncation  toward  zero). 


•  As  in  the  scalar  lattice  filter  and  WDF  cases,  sensitivity  of  coefficient  quan- 
tization  can  be  minimized  by  properly  scaling  the  network  to  deliver  “maximum 
power  transfer"  at  frequencies  where  low  sensitivity  is  required  [50j. 

•  A  perfectly  lossless  digital  realization  can  be  implemented  using  a  number 
system  presented  in  this  paper. 

•  The  desirable  structural  properties  are  derived  for  multi-input,  multi-output 
(MIMO)  transfer-function  matrices. 

The  derivation  of  the  WGF  is  made  exceedingly  simple  by  using  three  simple 
principles  of  wave  propagation  in  an  ideal  linear  medium.  To  our  knowledge,  these 
principles  have  not  been  invoiced  before  to  derive  digital  filters.  In  this  respect,  we 
feel  this  paper  has  significant  tutorial  value.  It  is  a  new  point  of  view. 

2.  Related  Prior  Work 

This  section  reviews  some  of  the  most  closely  related  work  on  digital  filter 
structures.  These  include  the  orthogonal- polynomial  filters  of  Szego,  the  “wave 
digital  filters"  of  Fettweis,  and  the  “orthogonal  filters"  of  Dewilde.  Naturally,  there 
are  many  more  related  lines  of  development,  in  view  of  the  first  law  of  signal 
processing.*  These  represent  only  the  major  recent  areas  closest  to  our  point  of 
view. 

2.1.  Wave  Digital  Filters 

The  wave  digital  filter  (WDF)  approach  of  Fettweis  (15,17,20,28]  comes  closest 
to  the  point  of  view  taken  in  this  paper.  Fettweis  obtains  a  similar  class  of  structures 
by  use  of  the  classical  notion  of  wave  variablet  [13j. 

*  The  first  law  of  signal  processing  is  “Everything  is  equivalent  to  everything  else." 


For  example,  if  v  and  »  denote  the  voltage  and  current  at  a  terminal  of  an 
N-port  network,  the  wave  variables  are  defined  by  r  =  v  +  Ri  and  y  s s  v  —  Ri, 
where  Aisan  arbitrary  "reference  impedance.”  These  wave  variables  are  logically 
equivalent  to  the  left-going  and  right-going  "pressure  traveling  waves”  considered 
in  this  paper,  and  R  plays  the  role  of  characteristic  impedance  in  the  associated 
transmission  line.  A  generalization  of  wave  variables  to  the  form  x  =  av+0i,  y  = 
TT«  +  d t  and  a  characterization  of  the  specialization  necessary  to  ensure  realizability 
is  given  in  [32|. 

Fettweis  describes  how  to  directly  model  resistors,  capacitors,  inductors,  trans¬ 
formers,  gyrators,  and  circulators  using  the  WDF  approach,  and  he  describes  the 
necessary  rules  for  connecting  ports  together  [17].  The  modeling  of  a  capacitor,  for 
example,  is  accomplished  by  scaling  the  reference  impedance  R  until  the  capacitor 
"reflectance”  is  exactly  a  unit-sample  delay.  (The  model  is  parametrized  in  fre¬ 
quency  so  that  the  wave  variables  are  really  phasors.)  An  inductor  also  maps  to  a 
unit-sample  delay  but  with  a  sign-change  relative  to  a  capacitor.  A  complete  circuit 
is  built  out  of  basic  elements  by  means  of  “adaptors”  (28|  which  play  the  role  of 
the  junctions  or  scattering  layers  described  in  this  paper;  the  adaptor  accomplishes 
interconnection  of  ports  at  different  reference  impedances. 

The  WDF  modeling  of  inductors  and  capacitors  is  limited  because  the  continuous¬ 
time  frequency  variable  is  mapped  to  the  discrete-time  frequency  variable  via  the 
bilinear  transform.  If  the  points  z  =  l  and  z  =*  —1  in  the  complex  plane  are 
identified  with  zero  and  infinite  continuous-time  frequencies,  respectively,  then  only 
one  more  mapping  frequency  (say  iii)  can  be  chosen.  Thus,  the  bilinear  transforma¬ 
tion  provides  exact  modeling  only  at  the  three  frequencies  0,  0,  and  oo. 

The  WDF  formulation  models  a  system  of  differential  equations  at  three  fre¬ 
quencies,  while  our  approach  exactly  models  wave  propagation  in  lossless  media 
having  spatially  discrete  changes  in  characteristic  impedance.  Consequently,  in  our 
formulation,  a  wave  variable  may  be  a  "voltage”  or  "current”  or  a  linear  combina¬ 
tion  of  the  two  without  incurring  realizability  problems  [32].  This  is  a  considerable 
conceptual  simplification  for  applications  to  physical  modeling. 


A  general  result  in  this  paper  is  that  overflow  oscillations  and  limit  cycles  can  be 
suppressed  in  a U  forms  of  scattering-type  filter  structures  simply  by  using  extended 
numerical  precision  in  each  scattering  section,  saving  quantization  (toward  zero)  for 
the  final  outgoing  waves.  The  basic  principles  involved  apparently  appeared  first  in 
the  WDF  context  (28,39). 


2.2.  Ladder  and  Lattice  Filters 


For  some  time  it  has  been  known  that  lattice  and  ladder  filtering  structures 
are  superior  to  the  so-called  direct  form  in  several  ways.  These  include  reduced 
sensitivity  to  coefficient  quantization,  less  dependency  of  round-off  noise  on  the  filter 
frequency  response,  ease  of  stability  checking,  reduced  probability  of  limit  cycles  or 
overflow  oscillations,  and  section-wise  orthogonality  in  the  linear  prediction  context. 
For  a  discussion  of  ladder  and  lattice  filters  in  adaptive  estimation,  see  [45], 

Lattice  structures  have  been  in  use  for  decades  in  directly  modeling  layered 
scattering  media.  The  mapping  of  underground  striations  in  rock  density,  for 
example,  is  a  basic  diagnostic  tool  in  oil  exploration.  The  interface  between  two 
subterranean  layers  of  rock  of  different  densities  produces  a  scattering  layer  because 
the  characteristic  impedance  of  the  medium  with  respect  to  sound  propagation 
changes  across  such  a  boundary. 

Another  example  of  the  use  of  lattice  structures  for  physical  modeling  is  the 
“acoustic  tube”  models  developed  for  speech  analysis  and  synthesis.  In  this  case, 
the  vocal  tract  is  modeled  as  a  cascade  of  coaxial  cylindrical  tubes  with  varying 
cross-sectional  areas  and  equal  length.  The  change  in  area  from  one  tube  section 
to  the  next  provides  a  change  in  the  characteristic  impedance  of  the  air  column  for 
sound  propagation,  and  so  a  series  of  equally  spaced  scattering  layers  is  obtained. 


Apparently,  the  filter  structures  developed  in  the  above  applications  are  only 
as  general  as  a  single  chain  of  scattering  layers  with  one  input  and  one  output,  and 
the  input  and  output  sections  are  terminated  in  a  non-extendable  way.  Little  if  any 


work  has  explored  branching  and  intersecting  chains  of  scattering  layers.  In  the  case 
of  speech,  the  use  of  a  separate  acoustic  tube  branching  off  from  the  vocal  tract  to 
mode!  the  nasal  tract  would  obviously  be  very  natural.  Apart  from  branching,  it  is 
not  possible  to  continue  the  structures  in  common  use  from  the  output  section  to 
a  larger  section.  This  is  because  the  typical  arrangement  is  to  assume  a  perfectly 
reflecting  termination  at  the  output.  Doing  this  allows  manipulation  of  the  delays 
in  the  scattering  network  to  place  them  more  conveniently  and  combine  in  pairs 
such  that  the  required  signal  sampling  rate  is  reduced  by  a  factor  of  2.  We  have 
found  that  the  cascade  scattering  chains,  which  dominate  the  recent  literature,  can 
be  immediately  extended  to  general  acyclic  trees  with  the  same  basic  properties. 

Our  formulation  is  more  general  than  even  the  acyclic-tree  extension  of  prevalent 
lattice  filters  in  that  arbitrary  networks  can  be  constructed.  Also,  there  does  not 
seem  to  be  an  existing  treatment  of  MIMO  systems  from  the  acoustic  waveguide 
point  of  view,  nor  the  generalization  which  allows  transmission  zeros  to  be  an  in¬ 
tegral  part  of  the  waveguide  (without  having  to  add  external  “taps"  for  forming  a 
linear  combination  of  the  contents  of  each  waveguide  section). 

A  particularly  important  antecedent  to  the  WGF  in  the  speech  processing 
literature  is  the  normalized  ladder  filter  (NLF)  developed  by  Gray  and  Marlcel 
[23,30,39].  Gray  considered  only  the  single-input,  single-output  (SISO)  all-pole 
case.  (Zeros  are  obtained  in  the  NLF  using  “taps,"  which  leads  outside  the  class  of 
structures  considered  here.)  Their  approach  was  based  on  orthonormal-polynomial 
expansion  [1,2,8]  which  is  closely  connected  with  linear  prediction  theory  [31].  They 
showed  the  following  to  be  true: 

•  The  NLF  is  optimal  in  the  sense  that  each  internal  node  has  unity  power 
gain.  This  means,  for  example,  that  the  response  to  a  unit  impulse  cannot 
overflow  anywhere  within  a  stable  NLF  filter.  Also,  if  the  input  signal  is 
white  noise  with  unit  variance,  the  variance  of  the  signal  at  each  internal 
node  is  exactly  unity  [30]. 

•  The  NLF  is  stable  in  the  case  of  time-varying  filter  parameters  [30]  as 
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long  as  the  "reflection  coefficients”  are  always  less  than  or  equal  to 
some  K<1  in  magnitude.  ([J:,(<)|  <  l  is  not  sufficient  for  bounded-input, 
bounded-output  (BIBO)  stability  unless  the  input  signal  energy  is  finite.)  It 
was  derived  incidentally  that  the  total  energy  entering  the  ladder  eventually 
"exits”  through  the  particular  delay  element  at  the  entrance  to  the  ladder. 

•  The  NLF  has  superior  roundoff  noise  properties,  especially  when  poles  are 
clustered  close  together  and/or  close  to  the  unit  circle  [30]. 

e  The  NLF  is  free  of  zero-input  overflow  oscillations  [39j. 

•  The  NLF  is  free  of  zeroinput  limit  cycles  [39]  in  magnitude-truncation 
arithmetic. 

The  NLF  is  obtainable  by  transformations  of  a  special  case  of  the  WGF  struc¬ 
tures  derived  here.  The  most  significant  difference  is  in  the  distribution  of  delay 
elements.  We  will  show  that  delay  distribution  in  the  standard  NLF  is  not  obtain¬ 
able  from  a  WGF  unless  the  waveguide  is  terminated  by  a  pure  reflection.  This 
means,  for  example,  that  an  NLF  cannot  be  connected  to  another  NLF  to  build 
a  larger  waveguide  system  with  finite  loading  from  one  stage  to  the  next.  Also, 
the  delay  distribution  chosen  for  the  NLF  is  such  that  creating  a  loop  with  NLF's 
yields  a  degenerate  (non-computable)  structure  because  a  delay-free  loop  appears. 
Another  limitation  of  the  NLF  is  that  the  concept  of  instantaneous  power  becomes 
artificial  for  individual  sections  (although  Gray  defines  a  non-physical  but  similar 
quantity  in  [39,  eq.  (2)[). 

A  disadvantage  of  the  NLF  is  that  it  requires  four  multiplications  per  pole  of  the 
filter  transfer  function.*  The  one-multiplier  lattice  filter,  on  the  other  hand  [31] 
requires  only  one  multiplication  per  pole.  Because  of  the  choice  of  delay  distribu¬ 
tion  in  standard  lattice-filter  theory,  only  the  NLF  has  been  shown  to  be  power¬ 
preserving  in  some  sense.  In  contrast,  we  will  show  that  even  our  counterpart 
to  the  one-multiplier  lattice  can  be  made  normalized  with  respect  to  time-varying 

*  As  a  side  result,  we  show  that  one  of  these  four  multiplications  can  be  eliminated. 


coefficients.  Conceptually,  this  is  achieved  by  compensating  the  amplitude  of  stored 
signal  samples.  The  resulting  normalized  one- multiplier  lattice  section  is  computa¬ 
tionally  less  expensive  than  the  NLF. 

The  reason  that  most  standard  ladder  and  lattice  structures  (all  but  the  NLF) 
cannot  be  power-normalized  in  the  time-varying  case  is  that  the  unnatural  distribu¬ 
tion  of  delays  adopted  makes  passivity  of  a  section  nontrivial  to  show.  This  paper 
describes  how  power- normalization,  perfect  energy  conservation,  and  complete  sup¬ 
pression  of  limit  cycles  and  overflow  oscillations  can  be  guaranteed  for  MIMO 
analogues  of  all  ladder  and  lattice  filter  structures,  with  extensions  to  branching 
structures  and  general  terminations. 

For  the  case  of  reflectively  terminated,  time-varying,  MIMO,  acyclic  trees, 

(which  specialize  to  ordinary  lattice/ladder  structures  in  the  SISO  single-branch 
case),  we  derive  efficient  equivalent  structures  in  which  the  delays  are  moved  and 
combined  to  yield  computational  savings  without  loss  of  the  desired  power-invariance/numerical 
properties. 

2.3.  Synthesis  and  Approximation 

The  synthesis  procedure  we  use  for  the  WGF  is  based  on  the  Schur  algorithm 
which  recursively  computes  a  solution  to  the  Nevanlinna-Pick  problem  [40,37,43]. 

The  Nevanlinna-Pick  problem  is  to  interpolate  a  rational  “Schur  function*”  through 
n  complex  values  at  n  points  in  the  closed  unit  disk  in  the  complex  plane.  The 
Schur  algorithm  has  also  been  called  the  “Nevanlinna  recursion  scheme”  (43).  In 
other  contexts,  a  special  case  of  the  the  Schur  algorithm,  which  computes  only 
all-pole  digital  filters,  has  been  called  the  “Durbin”  (8)  or  “Levinson”  [3]  algorithm 
(34,40,42,38,31).  The  complete  Schur  algorithm  constructs  a  cascade  WGF  realiza¬ 
tion  of  a  digital  filter  containing  both  poles  and  zeros. 

*  A  Schur  function  S(z)  is  defined  as  a  complex  function  analytic  and  of  modulus  not 
exceeding  unity  in  ]<(  <  1 


The  estimation  problem  has  been  addressed  by  DeWilde  [40,42].  In  this  cootext, 
the  Schur  algorithm  provides  an  ARMA  estimation  technique  in  which  the  pole 
estimates  are  optimal  in  the  mean  square  sense  for  the  given  fixed  zeros  which  are 
chosen  a  priori. 

3.  Traveling  Waves  and  Lossless  Scattering 

For  concreteness  of  discussion,  we  will  focus  on  pressure  and  flow  waves  in  a 
so-called  acoustic  tube.  We  could  just  as  easily  think  of  the  electric  and  magnetic 
components  of  light,  voltage  and  current  in  a  transmission  line,  or  force  and  trans¬ 
verse  velocity  on  a  vibrating  string.  An  analysis  of  the  acoustic  tube  is  discussed 
by  Market  and  Gray  [31}  and  Flanagan  [10]  in  the  context  of  vocal-tract  modeling. 
Further  details  on  the  acoustics  of  sound  in  tubes  can  be  found  in  Morse  [4|.  The 
term  “waveguide”  will  be  used  interchangeably  with  “acoustic  tube.” 

A  derivation  of  traveling  waves  from  the  basic  wave  equation  is  presented  in  the 
appendix.  The  result  is  that  in  a  cylindrical  acoustic  tube,  longitudinal*  pressure 
and  flow  waves  propagate  back  and  forth  with  speed  c.  Let  z  denote  distance 
along  the  tube  axis  and  let  t  denote  time  in  seconds.  Then  the  instantaneous 
pressure  P[z ,  <)  and  flow  £/(x,  t)  is  given  by  the  sum  of  the  left-going  and  right-going 
traveling-wave  components: 

P{z,t)=*P+(z,t)  +  P~(z,t)  (la) 

U(z,t)  =  U*(z,t)  +  U~(z,t)  (lb) 

*  We  assume  the  tube  radius  is  much  smaller  than  the  wavelength  of  sound  la  the  tube, 
so  that  pressure  and  flow  are  constant  over  any  cross-section  of  the  tube  normal  to  the 
axis.  In  other  words,  waves  do  not  propagate  up  and  down  but  only  left  and  right.  For 
more  details  on  the  assumptions  involved  in  acoustic  tube  models,  see  Flanagan  [10]. 


3.1.  Three  Fundamental  Constraints 


The  behavior  of  waves  traveling  unidirection  ally  in  a  lossless  medium  is  governed 
by  three  laws:  (1)  the  pressure  is  proportional  to  Sow,  (2)  the  pressure  is  a  continuous 
function  of  position,  and  (3)  the  flow  variable  (e.g.  mass  or  charge)  is  neither  created 
nor  destroyed  in  the  medium. 


Characteristic  Impedance 

An  ideal  linear  propagation  medium  is  completely  determined  by  its  charac- 
teristic  impedance^  Z(x,  t).  The  characteristic  impedance  is  defined  the  constant  of 
proportionality  between  pressure  and  flow  in  a  unidirectional  traveling  wave: 


P*  -  ZU+ 

(2a) 

P~  =  ~zu~ 

(2b) 

When  the  arguments  (r,  t)  are  omitted,  it  is  understood  that  all  quantities  are 
written  for  some  constant  time  l  and  position  z.  The  minus  sign  for  the  left-going 

'  For  an  acoustic  tube,  the  characteristic  impedance  is  given  b jr  Z  ■  JPotTp/S  »■> 
pc/S,  where  p  is  the  density  (mass  per  unit  volume)  of  air  in  the  tube,  e  is  the  speed  of 
propagation,  Pc  is  ambient  pressure,  fe  is  the  ratio  of  the  specific  heat  of  air  at  constant 
pressure  to  that  at  constant  volume,  and  5  is  the  cross-sectional  area  of  the  tube.  In  a 
vibrating  string,  Z  «■  JTp  —  pc,  where  p  is  string  density  (mass  per  unit  length)  and  T  is 
the  tension  of  the  string.  In  an  electric  transmission  line,  Z  ■■  \jLjC  “  Lc  where  L  and 
C  are  the  inductance  and  capacitance,  respectively,  per  unit  length  along  the  transmission 
line.  In  free  space,  Z  ™  \J po/to  “  poc,  where  p«  and  ««  are  the  permeability  and 
permittivity,  respectively,  of  free  space. 
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wave  P~  accounts  for  the  fact  that  flows  in  opposite  directions  subtract  while 
pressure  waves  passing  through  each  other  add. 

We  will  consider  initially  a  more  general  situation  in  which  Z  =  Z{d)  is  a 
q  by  q  complex  matrix  function  of  the  complex  variable  (or  unit-delay  operator)  d. 
For  stability  of  propagation  in  the  waveguide,  we  require  that  Z(d)  be  analytic  for 
\d\  <  l.  The  results  also  extend  to  the  case  of  vector  &  ==  [d\, . ..  ,ije],  &ut  we 
will  treat  only  one  complex  argument  d  for  notational  simplicity.  The  pressure  and 
flow  variables  are  q  by  m  matrix  complex  functions  of  d.  However,  keep  in  mind 
that  the  physical  analogy  we  are  pursuing  is  for  the  case  of  real  scalar  Z,  P,  and 
U. 

For  lossless  propagation  in  the  scalar  case,  the  characteristic  impedance  Z 
must  be  real.  In  the  matrix-delay-operator  case,  lossless  propagation  will  now  be 
characterized  by  the  requirement  that  Z  be  para-Hermitian,  i.e., 

Z.(d)  *  Z{d)  (3) 

where 

Z.(d)AT( 1/2?  (4) 

denotes  the  para-Hermitian  conjugate  of  Z(d)  (I3,40j,  (-)r  denotes  transposition, 
and  denotes  complex  conjugation.  For  d  »  e*9,  Z,{c*9)  coincides  with  the 
Hermitian  transpose  of  Z[e}9).  The  para-Hermitian  conjugate  is  the  unique  analytic 
continuation  (when  it  exists)  of  the  Hermitian  transpose  Z«(e**)  =  Z(eJ9)  from 
the  unit  circle  into  the  complex  plane.  Thus,  a  lossless  medium  in  our  framework 
is  defined  as  one  in  which  the  characteristic  impedance  is  para-Hermitian.  The 
extension  to  vector  i  is  obtained  by  regarding  Z(i)  as  K  functions  of  scalar  complex 
variables  df.  Note  that  in  the  scalar  case,  Z  para-Hermitian  impties  Z  -=Z  which 
implies  Z  is  real.  Henceforth,  we  assume  Z  denotes  a  para-Hermitian  characteristic 
impedance.  For  non-para-Hermitian  Z,  (2)  should  be  modified  to  read  P  = 
— Z»U~  (13],  and  a  passive  medium  is  one  in  which  Z  +  Z.  is  positive  semi-definite. 


It  is  worthwhile  to  interpret  the  various  levels  of  extension  we  are  consider¬ 
ing  for  the  characteristic  impedance  Z.  When  Z  is  real  and  scalar,  we  obtain  ex¬ 
actly  the  ideal  behavior  of  one  dimensional  traveling  waves  in  a  lossless  medium. 
Extending  to  q  by  q  matrix  characteristic  impedances  facilitates  development  of 
multi-input,  multi-output  (MIMO)  systems  which  have  the  desired  numerical  and 
power-invariance  properties.  The  extension  to  analytic  matrix  functions  of  a  com¬ 
plex  variable  provides  a  generalized  scattering  medium  whose  reflectance  and  tran»* 
mission  coefficients  are  themselves  rational  transfer  function  matrices.  This  provides 
for  nesting  of  the  WGF  structures.  The  complex  argument  d  of  the  characteristic 
impedance  is  interpreted  as  a  unit-delay  operator,  and  the  meaning  of  the  charac¬ 
teristic  impedance  is  attached  to  its  Laurent  series  expansion  with  respect  to  the 
unit  circle  in  the  d- plane.  Additional  complex  variables  d,  in  the  arguments  to  the 
characteristic  impedance  allow  the  generalized  scattering  layer  to  perform  filtering 
in  several  domains  such  as  time  and  space.  Since  the  characteristic  impedance  is 
assumed  stable  and  para-Hermitian,  all  delay-operator  impedance  matrices  must 
be  nonrecursive  and  zero-phase.  Therefore,  computability,  stability,  and  nonlinear 
oscillation  problems  do  not  arise  in  the  case  of  multiple  domains. 

Pressure  Coatiauity  a  ad  Medium  Co  as  emtio  a 

We  will  be  interested  in  the  situation  wherein  the  characteristic  impedance 
changes  abruptly  from  one  value  to  another,  say  from  Z\  to  Z3.  The  impedance 
discontinuity  can  be  a  sudden  change  along  r  in  the  acoustic  tube,  or  it  can  be 
a  change  introduced  at  some  time  t  (as  needed  for  time- varying  filters).  First  we 
consider  changes  with  respect  to  r.  Given  the  traveling  waves  impinging  on  the 
junction  between  Z\  and  Z2,  we  seek  formulas  for  the  traveling  waves  leaving  the 
junction.  To  solve  this  problem,  we  need  two  laws  in  addition  to  (2)  for  an  ideal 
wave  medium: 

1)  Pressure  cannot  change  instantaneously  across  the  junction  (5a) 

2)  The  sum  of  flows  meeting  at  the  junction  is  zero  (5b) 


WJ  1  1  \m  I  WJ  Af'  ■  H.I  I  1.  1  KJ  mi  gi»  ( 
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These  constraints  are  often  called  *Kirchoff’s  node  equations”  in  the  context  of 
c;-cuit  theory.  For  changes  in  characteristic  impedance  with  respect  to  time,  (S) 
is  not  applicable.  Time-varying  characteristic  impedances  will  be  implemented  as 
waveguide  transformers,  and  will  be  used  to  obtain  power-invariant  lossless  digital 
filters  in  the  time-varying  case. 

The  continuity  and  conservation  constraints  (5)  together  with  the  characteristic 
impedance  constraints  (2)  determine  what  pressure  and  flow  waves  emerge  from  a 
junction  between  waveguide  media  of  differing  characteristic  impedance,  given  the 
incoming  waves. 

Consider  the  case  of  N  waveguides  meeting  at  a  common  junction.  Kirchoff's 
laws  state  that  there  can  be  only  one  resultant  pressure  Pj  at  the  junction,  and  the 
sum  of  flows  entering  and  leaving  the  junction  must  total  to  zero.  Thus,  we  have 
the  constraints 


Pi  =  P2  -  • 

Ui  +  Ui  + 


=  Pn 

•  +  UN 


Pj 

0 


(6a) 

(6b) 


where 


p*  -  P*  +  Pj  p!  =  ZiU- 

Ui  *  uj  +  U~  p-  m  -ZiU+i 

Characteristic  impedance  of  the  ith  waveguide  (q  by  q) 


(7) 


Zi 

r,-  »  Z~l  »  Characteristic  admittance  of  the  ith  waveguide  (g  by  7) 

Pi  *  Incoming  pressure  wave  along  the  ith  waveguide  (q  by  m ) 

U i  *  Incoming  flow  wave  along  the  ith  waveguide  ( q  by  m) 

P{  »  Outgoing  pressure  wave  along  the  ith  waveguide  (q  by  m) 

Uj  =  Outgoing  flow  wave  along  the  ith  waveguide  (q  by  m) 

Pi  *  Instantaneous  pressure  wave  in  ith  waveguide  just  outside  junction  (q  by 
Ut  Instantaneous  pressure  wave  in  ith  waveguide  just  outside  junction  (q  by 
Pj  a  Resultant  pressure  in  the  junction  (q  by  m) 


(S) 


m) 

m) 

(9) 
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3.2.  Reflection  at  a  Junction 


Given  a  set  of  incoming;  traveling  pressure  waves,  P* ,  i  =  the 

constraints  (2,6)  determine  the  outgoing  waves  as  follows.  As  before,  Z ,•  and  hence 
T,-  »  Z~l  are  para-Hermitian  positive  definite.  Substituting  (lb)  and  (2)  into  (6b) 
yields  the  resultant  junction  pressure: 


(N  \_t  .V 

Er-j  E  r.-e; 


r j  A  Er'  Zj&n1  <u> 

•—I  1—1 

define  the  junction  admittance,  junction  impedance,  and  junction  flow,  respectively. 
(In  the  extension  to  non- para-Hermitian  Uj  becomes  U*j  =  ]T(r,-  +  I\«)P*.) 
Relation  (10)  then  reads  Pj  =  ZjUj,  or. 


Junction  Pressure  =  Junction  Impedance  x  Junction  Flow 


Since  I\P*  =  U* ,  we  have  Uj  =  2££_,  U*  &  2U*  =»  |C'yl  =  ICr+|+|C:”l, 
where  |  •  |  denotes  elementwise  complex  modulus.  That  is,  the  junction  flow  can 
be  interpreted  as  the  magnitude  sum  of  the  incoming  and  outgoing  current  waves. 

Now,  given  incoming  traveling  waves  P* ,  U*  and  the  characteristic  impedance 
Z(  of  each  branch  terminating  at  the  junction,  we  easily  find  the  outgoing  waves 
Pj,  U~  to  be 

P~  =Pj-  P*  (12a) 

uj  =  r  tPj  (12b) 

Equations  (12)  specify  the  scattering  at  the  junction  of  N  intersecting  “wave  guides,” 
given  the  incoming  waves  P+  (or  U*)  and  the  branch  characteristic  impedances  Z,, 
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In  view  of  relations  (2),  we  can  consider  only  left-going  and  right-going  pres¬ 
sure  waves,  since  the  flow  waves  can  be  readily  computed  from  the  characteristic 
impedance  of  the  propagation  medium.  At  this  point,  we  could  instead  choose  wave 
variables  of  the  form  r,-  =  P,  +  Z,C/,-  and  y;  =  P,  —  Z,C/,-  and  proceed  along  the 
lines  of  classical  wave  filtering  (13].  However,  such  a  path  is  less  fundamental  in  the 
present  development  because  we  are  considering  only  discrete-time  filters. 

We  have  treated  only  a  parallel  junction  of  waveguides.  A  dual  set  of  equations 
is  obtained  by  considering  a  aeries  junction.  However,  pressure  waves  intersecting 
in  a  parallel  junction  are  equivalent  to  flow  waves  intersecting  at  a  series  junction. 
When  using  flow  waves  as  the  primary  variables,  (12b)  can  be  written 

V~  =>  V'  -  r,Pj 

Pj  =  iZj  £  u~ 

lol 

The  series  pressure  junction  is  obtained  by  talcing  the  dual  of  (13).  That  is,  replace 
£/,  by  P,  and  T,  by  Z,  to  obtain 


p r = pf  -  ZiVj 

(Hi) 

Vj  =  2  n  £  p; 

(14b) 

n  -  zrr' 

(14c) 

<v 

Z"j  =  E  z, 

1—1 

(14d) 

The  junction  impedance  for  a  series  junction  is  the  sum  of  the  branch  impedances, 
while  for  a  parallel  junction,  it  is  the  parallel  combination  of  the  branch  impedances 
(inverse  of  the  sum  of  admittances). 

Equation  (12a)  is  a  computationally  efficient  way  to  implement  an  Alport 
scattering  junction.  In  the  case  iV  =  2,  the  well-known  one-multiplier  lattice  filter 
section  (minus  its  unit  delay)  is  obtained  immediately  from  (12a).  More  generally. 


(13a) 
( 13b  j 
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40  iV-way  intersection  requires  N  multiplies  and  N  —  1  additions  to  obtain  Pj, 
and  one  addition  for  each  outgoing  wave,  for  a  total  of  N  multiplies  and  2 N  —  1 
additions.  The  dual  junction  (14)  also  requires  N  multiplies  and  2 N  -  1  additions. 
In  the  next  section,  a  method  for  trading  one  multiplier  for  another  IV—  1  additions 
[28]  is  described. 

3.3.  o-parameters 


One  parametrization  of  all  passive  N-junctions  is  the  set  of  N  branch  im¬ 
pedances  with  positive-definite  para-Hermitian  parts  (cf.  §5.l).  This  section  describes 
another  parametrization,  analogous  to  that  used  in  the  WDF  context  [28]. 

One  parametrization  of  all  passive  jV-junctions  is  the  set  of  i\  branch  im¬ 
pedances  with  positive-definite  para-Hermitian  parts.  This  section  describes  another 
parametrization,  analogous  to  that  used  in  the  WDF  context  [28]. 

Define 

or,  »  2Zyr,  (15) 

which  is  twice  the  junction  impedance  times  the  ith  branch  admittance.  (In  the 
non-para-Hermitian  case,  or,-  =  Zy(r,-  +  r ,-.).)  Then  the  junction  pressure  can  be 
written  as  a  linear  combination  of  the  incoming  pressure  waves  in  terms  of  the  a, 

as 

N 

Pj  =  tTq<p>  fI6) 

i—i 

Sine.  Eli,  r.  4r/t 

N 

£ar,-  =  2/,  (17) 

•—1 


where  Iq  is  the  q  by  q  identity  matrix. 


In  matrix  form,  (12a)  can  be  written 


or 


where 


* 

’  dl  —  [q  ...  <*/V 

at  a2  ~  ... 

/>* 

m 

•  •  • 

.  •  • 

• 

p* 

at  a2  •••  a/v  — 

• 

£ 


(18) 


(19) 


(20) 


The  matrix  E  is  called  the  scattering  matrixot  the  junction.  Since  P~  =  (a,— lq)P+ 
when  Pj  =  0  for  all  j  i,  we  define  the  reflection  coefficient  at  the  «th  port  by 


(21) 


Equations  (12a, 16, 17)  combine  to  give 

(V 

p7  =  P*  +  E  Q;<p/  -  p>  )  (22) 

;  — i 

iri' 

Thus,  ay  can  be  interpreted  as  the  fraction  of  the  pressure  differential  between 
branches  j  and  «  which  is  reflected  back  along  the  ith  branch  with  P* ,  for  any  ». 
Use  of  this  expression  saves  one  matrix  multiply  but  entails  3N—2  matrix  additions. 
If  one  multiply  is  worth  Af  —  1  additions  or  more,  then  (22)  is  less  expensive  to 
implement  than  (12a). 
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3.4.  Pure  Delays 


Up  to  now  we  have  been  concerned  only  with  the  scattering  of  traveling  waves 
at  a  single  impedance-dbcontinuity  junction.  We  now  allow  for  many  such  junctions 
interconnected  by  waveguides  for  lossless,  reflectionless  propagation.  Physically,  an 
interconnection  between  junctions  is  a  length  of  material  at  a  single  characteristic 
impedance.  It  is  implemented  digitally  using  bi-directional  delay  lines. 

Consider  the  interconnection  of  two  iV-port  junctions.  Between  the  two  junc¬ 
tions  is  a  section  of  pure  waveguide  which  is  a  lossless  medium  having  characteristic 
impedance  Z jj.  Let  c  denote  the  speed  of  propagation  in  this  waveguide  section, 
and  suppose  the  distance  between  the  junctions  is  L.  Then  the  propagation  time 
from  one  junction  to  the  other  and  back  is  7*.  —  2 L/c.  Consider  a  pressure  wave 
impulse  P*(x  —  ct)  =  6(x  —  ct)  traveling  from  junction  1  to  junction  2  starting 
at  time  t  =  0.  At  time  T,/ 2  it  reaches  junctions,  and  a  reflection  P  (x  +  cf)  = 
/»l36[x  +  c(t  — T,)!  starts  out  to  the  left  from  junction  2,  heading  back  to  junction  1. 
A  fragment  of  the  pulse  is  also  sent  out  along  all  waveguides  connected  to  junction 
2,  according  to  the  relative  branch  impedances.  At  time  t  =  T,,  the  reflected  pulse 
reaches  junction  l,  and  scatters  away  again. 

A  section  of  waveguide  joining  two  junctions  by  a  propagation  delay  T,[2  is 
called  a  unit  delay.  If  the  speed  of  propagation  is  everywhere  the  same,  then  all 
unit  delays  are  of  the  same  length  L  =  cT,/ 2. 

An  important  observation  to  make  at  this  point  is  that  the  impulse  response 
at  a  particular  junction  of  a  network  of  junctions  interconnected  by  waveguides  of 
length  cTt/2  b  nonzero  only  at  times  which  are  integer  multiples  of  Te.  Thus,  such 
a  network  b  a  simulation  of  a  digital  linear  system  with  sampling  rate  Ft  =  1  /Tt. 

Consider  a  linear  chain  of  junctions.  The  input  b  defined  as  the  pressure  enter¬ 
ing  at  the  far  left  (junction  0),  and  the  output  b  defined  as  the  pressure  emerging 
from  the  far  right  (junction  M).  Again,  thb  structure  b  an  exact  simulation  of  a 
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digital  system  with  sampling  rate  F «  =s  1/T*  provided  M  is  even.  If  M  is  odd,  we 
obtain  samples  separated  by  T,  seconds  at  the  output,  but  there  is  a  time  shift  of 
T«/2  relative  to  the  sample  instants  at  the  input. 

Now  consider  the  more  elaborate  interconnection  of  junctions  in  an  arbitrary 
network.  Unless  every  path  from  the  input  to  the  output  crosses  an  even  number 
of  unit  delays,  the  impulse  response  of  the  network  will  be  nonzero  at  multiples  of 
Tt/ 2,  yet  the  minimum  delay  in  any  feedback  loop  is  Tt  seconds.  Thus,  it  is  not 
possible  to  simulate  a  general  digital  system  having  sampling  rate  2/7*,  even  though 
such  a  sampling  rate  is  required  to  compute  the  response. 

Two  cases  arise: 

(1)  Half-rate  structures  which  require  an  even  number  of  branches  on  every 

path  from  the  input  to  the  output,  plus  possibly  a  single  odd  section 
which  causes  a  half-sample  shift  of  the  output  relative  to  the  input.  In 
the  latter  case,  we  do  not  allow  the  resulting  structure  to  be  placed  in  a 
feedback  loop.  These  restrictions  leave  us  with  a  general  class  of  linear 
digital  filter  structures. 

(2)  Full-rate  structures  in  which  the  most  general  transfer  function  between 

two  junctions  is  of  the  form 


//(*)  - 


6t  Z~l  +  ••  •  4-  bnk  Z~"> 
1  +  0  aj  z~ 2  +  •  •  •  +  a„. 


in  the  scalar  case,  where  na  >  2  and  nj  >  1  are  integers.  In  the  full-rate 
case  we  must  settle  for  the  class  of  rational  filters  in  which  at  =  0.  The 
corresponding  digital  filter  design  methods  must  support  this  constraint. 
Therefore,  techniques  based  on  linear  programming  seem  natural  in  this 
context  (38|.  In  the  estimation  context,  the  methods  described  in  [49]  can 
be  adapted  to  this  problem. 
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4*  Cueade  Waveguides 


We  oow  specialize  discussion  to  the  case  of  cascade  waveguide  sections.  The 
junction  between  two  waveguides  of  differing  characteristic  impedance  will  create 
a  scattering  junction.  The  stretch  of  pure  waveguide  material  between  scattering 
junctions  will  provide  delay  lines  for  the  propagation  signal.  From  these  structures 
all  digital  filters  can  be  built  in  such  a  way  that  they  behave  nicely  with  respect 
to  time-varying  parameters  and  numerical  roundoff/overfiow.  Note  that  almost  all 
special  properties  in  the  cascade  case  carry  over  to  arbitrary  acyclic  trees. 


4.1.  The  Two-Port  Junction 

If  there  are  only  two  waveguides  meeting  at  a  junction,  we  obtain  the  classical 
“scattering  theory"  in  which  an  incoming  wave  is  split  into  a  “reflected”  and 
“transmitted"  part.  From  (15),  the  o-parameters  are 

c,  =  2(zr'  +  ZVY'Z-'  -  2(r,  +r,r‘r,  =  +  z, r,) 

(04) 

<,,  =  2(27' +  0  'zrl=2(r,  +  r,rIr,-j(/,  +  z2r1) 

From  (22),  the  reflected  pressure  waves  are 

P~  =  a2[p;  -  />r)  +  P*  -  (at  -  I<t)P*  +  ot2P n 

(°5) 

p 2  +  +(a2-/,)/>2 

If  P2  *  0,  then  the  incidence  of  P ^  produces  a  reflected  wave  P,  =  (a^  —  Iq)P*  ■ 
Thus,  we  define  the  following  reflection  coefficients: 

pi  AQl  -/,  =  (z2  -  zx)(z2  +  z,rl  =  (rt  +  r2r‘(rt  - r2) 

A  ,  (*s) 

P7  —  ®2  —  Iq  =  — Pi 


41 


It  is  now  apparent  that  if  the  reflection  coefficient  at  port  1  is  pi  —  p,  then  at 
port  2  the  reflection  coefficient  is  -p.  Another  point  of  view  is  that  inverting  the 
impedance-step  ratio  from  Zi  :Z^  to  Zq:Z\  merely  changes  the  sign  of  the  reflection 
coefficient.  It  b  easy  to  show  that  for  scalar  Z, ,  exchanging  pressure  waves  for  flow 
waves  also  toggles  the  sign  of  all  reflection  coefficients  in  the  network.  Thus,  left  b 
the  dual  of  right  as  pressure  b  the  dual  of  flow(and  parallel  b  the  dual  of  series). 

In  the  matrix-impedance  case  (q  >  I),  however,  replacement  of  pressure  by 
flow  changes  p  to  — (T i Zq  —  IqH^iZq  +  Iq)~l  which  equab  —  p  only  if  Z\  commutes 
with  Zj.  Two  Hermitian  matrices  commute  if  they  have  the  same  eigenvectors,  i.e., 
their  “principal  axes  of  dilation”  are  aligned.  There  exbts  a  unitary  transformation 
of  any  Hermitian  matrix  which  commutes  with  any  other  Hermitian  matrix;  that  b, 
a  Hermitian  matrix  can  be  “rotated"  until  it  commutes  with  any  other  Hermitian 
matrix.  For  waveguides  with  impedance  matrices  so  aligned  (all  Z,-  have  the  same 
eigenvectors  on  the  unit  circle),  junction  reversal  b  equivalent  to  wave-variable 
exchange;  either  causes  a  sign  reversal  in  all  reflection  coefficients. 


4.2.  The  Scattering  Matrix 

In  block-matrix  notation,  the  junction  output  is  given  by 


o.  -u  °»  livr 

o,2  a2-Iq  p* 


Pi  Iq  +  Pi 
Iq  -  Pi  -Pi 


A  (27) 


or 

Z  =  (2S) 

This  b  called  the  scattering  formulation,  and  E  b  called  the  scattering  matrix.  It 
b  a  special  case  of  the  jV-junction  scattering  matrix  defined  in  (18). 


4.3.  The  Chain-Scattering  Matrix 


In  the  two-port  case  only,  we  can  also  define  the  ehain-acattering  matrix  via 


KJ 


=  e 


Px 

p * 
r2 


(29) 


where 


0  A 


(30) 
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While  the  scattering  matrix  computes  outgoing  waves  from  incoming  waves,  the 
chain-scattering  matrix  computes  the  left-going  and  right-going  waves  in  section  l 
given  the  left-going  and  right-going  waves  in  section  2.  The  relation  between  the 
scattering  matrix  E  and  the  chain-scattering  matrix  6  is  given  by 


©11  —  -Ji  ©12  “  “S21  £22 

©21  —  Eu-2il  ©22  =  E  —  WSa 


(31) 


and 


Et|  =  ©21  ©1  /  E  —  ©22  —  ©2i©i/© 


E-i  =  ©r,1 
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(32) 


4.4.  One-Multiplier  Forms 


Equations  (25)  and  (26)  combine  to  give 

P~  =  P*  +  AP 
Pi  =  Pi  +  AP 


(33) 


Thus,  only  one  matrix  multiplication  is  necessary  to  compute  the  reflected  waves 
from  the  incoming  waves.  In  the  scalar  case,  this  reduces  to  the  so-called  one- 
multiplier  lattice  section  [31]  (minus  the  unit-sample  delay  ordinarily  associated 
with  each  section).  It  is  well-known  that  any  rational  digital  filter  can  be  built 
using  one-multiplier  lattice  sections  (31j.  In  fixed-point  implementations,  the  only 
source  of  error  would  typically  be  in  the  computation  of  A P. 

To  ensure  the  absence  of  limit  cycles  and  overflow  oscillations,  the  additions 
in  (33)  must  be  performed  before  rounding,  and  the  final  rounding  to  obtain  P1 
and  Pj  must  be  norm  reducing.  (In  the  scalar  case,  magnitude  truncation  is 
sufficient.)  The  added  expense  for  postponing  round-off  until  the  final  outgoing 
waves  are  computed  is  typically  negligible,  requiring  only  logic  to  determine  the 
desired  direction  of  truncation  from  the  signs  of  P*  and  AP,  and  the  low-order 
product  of  AP.  In  other  words,  the  double  precision  computations  required  are 
only  conceptual  because  the  low-order  half  of  the  incoming  waves  is  zero. 

Another  one-multiplier  form  is  obtained  by  organizing  (25)  as 


p:=p; +o(Pj  -p*) 

p;  =  pr  +  (pr-p;>  =  py-p; 


(35) 


where  a  &  at.  As  in  the  previous  case,  only  one  multiply  and  three  adds  are 
required  per  section.  An  advantage  to  (35)  is  that  the  computation  of  P2  is  noiseless 
(assuming  roundoff  only  after  multiplication).  This  can  be  used  to  simplify  hardware 
for  suppressing  limit  cycles. 

In  the  scalar  case,  the  single  section  parameter  p  of  (35)  must  lie  between  —  1 
and  1,  while  in  (34),  the  single  section  parameter  a  must  lie  between  0  and  2. 
Otherwise,  the  junction  is  not  passive.  The  practical  implication  of  non-passive 
junctions  is  potential  filter  instability  in  the  presence  of  feedback. 
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6.  Signal  Power  In  Lossless  Waveguides 

This  section  presents  basic  definitions  of  signal  energy  and  power  in  a  waveguide. 
These  concepts  provide  the  necessary  handles  on  filter  stability,  passivity  with 
respect  to  numerical  round-off  and  overflow,  and  energy  modulation  due  to  time- 
varying  parameters. 

5.1.  Instantaneous  Propagating  Power 

The  instantaneous  power  in  a  waveguide  containing  instantaneous  pressure  P 
and  flow  U  is  defined  as  the  m  by  m  product  of  pressure  and  flow: 

P  =  PmU  (36) 

The  total  instantaneous  power  is  defined  as  the  trace  of  the  instantaneous  power. 

Pp  ATr{P}  =Tr{P.f/}  (3~) 

The  total  instantaneous  power  is  a  complex  scalar  measure  of  power  flow.  It  can 
be  interpreted  in  a  manner  similar  to  the  complex  power  in  scalar  transmission¬ 
line  theory  in  which  sinusoidal  phasors  are  propagated  in  either  direction.  The 
instantaneous  power  can  be  expanded  into  four  terms  as  follows: 

P  =  P.U  =  ( P *  +  P~)(U *  +(/")  =  (U*  -  U~)Z.(U *  +  U~) 

=  P*U*  +  P~U~  +  P*U~  +  P~M:* 

.  -  ♦  (3*1 

=  pyp*  -  p~vp~  +  p, tp"  -  p.  vp 

=  U*ZU*  -  U~  ZU~  -  U* ZU~  +  1 7' zu* 


The  right-going  and  left-going  power  are  defined,  respectively,  by 

p*  *  p+u+  *  u*.zv*  =  P^rp^ 

P~  *  P~u~  =  -U~ZU~  =  -P~VP~ 


(30) 
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Since  Z  is  para-Hermitian,  P+  and  P  are  Hermitiaa  forms,  and  can  be  expressed 

as 

m 

P+  = 

/  i  i  “i* 

T  (40) 

>p"  =  V  X-u-u- 

/  _^«  I  "“I  •• 

I— l 

where  Of  is  the  tth  eigenvector  of  P+,  and  \*  is  its  »th  (real)  eigenvalue.  The 
m*vectors  U*  can  be  chosen  orthonormal.  Similar  remarks  apply  to  the  eigenvalues 
and  eigenvectors  of  P~ .  It  can  be  shown  that  the  waveguide  is  passive  if  and  only 
if  X*,X^  >  0.  Consequently,  we  will  assume  in  the  sequel  that 

xt  >0,  X-  >  0,  «'  =  1,2 . m  (41) 

This  implies  P*  and  P  are  positive  definite  Hermitian  matrices.  The  orthonormal 
vectors  and  U~  (which  are  vector  analytic  functions  of  <f),  indicate  “directions' 
along  which  power  flows  in  the  m«duneasional  manifold  determined  by  U,  (or  P,) 
and  Z. 

In  the  non-para-Hermitian  case,  the  medium  is  passive  iff 

Z[d)  +  Z.{<*)  >  0,  V|d|  <  l  (42) 

That  is,  the  para*Hermitian  part  of  the  characteristic  impedance  of  a  passive 
medium  must  be  positive  definite  in  the  unit  circle.  It  can  be  shown  using  the 
maximum  modulus  theorem  (5]  that  (42)  holds  if  Z  +  Z.  >  0  for  |d|  *  l. 

We  define  the  cross  power  by 

P*  =  U+ZU~  (43) 

The  instantaneous  power  (38)  can  now  be  written  as 

/>=(/>♦-/>-)  +  (/>*-/>*)  (44) 
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which  we  interpret  as  &  sum  of  &  net  traveling-power  term  P*  —  P  plus  the  skew- 
pararHermitian  part  of  the  cross-power.  In  the  real  scalar  case,  P  =  P.  and  the 
cross  power  is  zero. 

Since  the  eigenvalues  of  a  Hermitian  matrix  are  purely  real,  we  define  the 
difference  between  the  right-going  and  left-going  power  P*  —  P  as  the  active  power. 
Similarly,  since  the  eigenvalues  of  a  skew-Hermitian  matrix  are  purely  imaginary, 
we  define  the  skew-para-Hermitian  part  of  the  cross-power  P  —  P,  as  the  reactive 
power.  These  definitions  parallel  those  of  scalar  transmission- line  theory.  The  total 
power  in  each  case  is  defined  by  the  corresponding  trace. 

5.2.  Power  at  a  Junction 

For  the  A’-way  waveguide  junction,  the  constraints  (2,6)  yield 
N  N  ,V 

pjAZ  pm  -  E  pm.  -  pj  X>=°  («) 

i—l  »—l  i— t 

Thus,  the  N- wav  junction  is  lossless;  no  net  power,  active  or  reactive,  flows  into  or 
away  from  the  junction. 

5.3.  Quantization  Effects 

While  the  ideal  waveguide  junction  is  lossless,  finite  wordlength  effects,  can 
make  exactly  lossless  networks  unrealizable.  In  fixed-point  arithmetic,  the  product 
of  two  numbers  requires  more  bits  (in  general)  for  exact  representation  than  either  of 
the  multiplicands.  If  there  is  a  feedback  loop  around  a  junction,  the  number  of  bits 
needed  to  represent  exactly  a  circulating  signal  grows  without  bound.  Therefore, 
some  sort  of  round-off  rule  must  be  included  in  a  finite-precision  implementation  of  a 
WGF.  The  guaranteed  absence  of  limit  cycles  and  overflow  oscillations  is  tantamount 
to  ensuring  that  all  finite-wordlength  effects  result  in  power  absorption  at  each 
junction,  and  never  power  creation. 


The  sum  of  power  flows  entering  s  junction  b  given  by 


Pj-E.IL  A  £  />,.(/, 

»— l 

4  z(f>:.+r:.)(u:+v;) 

i—  I 

=  E(/’r.+p-)r.(p;-/>r) 

•—t 


where  T,  =  Z~l  is  the  characteristic  admittance  of  the  *th  waveguide. 

N 

i—i 

sV  - 

II*  llc  A  £*.-.r.-*.--,/(iflc 

•—  t 

Then  by  ( JO), 


+£-,£*-r)c 

lH£1c  +  (£''£V(£*'r)c 


Define 


The  junction  is  passive  if 


Since  T,-  >  0,  the  quadratic  form  x.r,*x  forms  an  elliptic  norm  for  r.  By  the 
norm  equivalence  theorem,  condition  (49)  can  be  written 

<5°> 


where  the  norm  is  arbitrary  (13].  Let  P{  =  Pj  —  e,-  denote  the  quantized  value 
of  P~.  Then  a  sufficient  condition  for  the  absence  of  limit  cycles  and  overflow 
oscillations  in  an  jV-port  junction  is 


<  \P 


(51) 


Since  the  norm  is  arbitrary,  two  convenient  choices  are  the  Ll  norm  (maximum 
abolute  column  sum)  and  the  L°°  norm  (maximum  abolute  row  sum)  [IS] .  Alternately, 
a  sufficient  condition  for  the  absence  of  overflow  oscillations  and  limit  cycles  in 
networks  built  from  N -port  waveguide  junctions  is.  that  magnitude  truncation  be 
used  on  each  element  of  the  final  q  by  m  outgoing  wave  variables  P . 


5.4.  The  Normalized  Ladder  Section 


We  can  normalize  the  pressure  and  flow  variables  by  the  Hermitian  square  root 
of  the  characteristic  impedance  to  obtain  propagation  waves  in  units  of  root  power: 

K  ±z’l'Pi  r~  ±*7*p? 

(52) 

0*  &z}u*  U~  AZ?U~ 

where 

A  =  A  (S3) 

is  the  unique  para-Hermitian  square  root  of  Z,-.  The  para-Hermitian  square  root  of 
Zi  is  defined  as  the  analytic  continuation  of  the  Hermitian  square  root  of  Zjic*0). 
Uniqueness  is  inherited  from  uniqueness  of  both  the  Hermitian  square  root  and  the 
process  of  analytic  continuation. 
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In  the  non-para-Hermitian  case,  we  normalize  the  traveling  waves  by 


P*  A  R~  •  P*  P~  A  R~  * Pj 

u*  Ar\u*  U~  A  r)u~ 


(54) 


where 

R.  k  ;(Z.  +  2,.) 

(55) 

1.1  '  ' 

Rl  =  R? 

i  . 

That  is,  R f  is  the  para-Hermitian  square  root  of  the  para-Hermitian  part  R,  of  the 
ith  branch  impedance  Z ,  (13]. 


By  restricting  all  waveguides  to  normalized  waves,  we  obtain  the  WGF  generaliza¬ 
tion  of  the  normalized  ladder  form  (NLF).  As  a  result  of  this  normalization,  the 
stored  power  in  each  WGF  sectioo  is  unaffected  by  time  variation  of  the  characteris¬ 
tic  impedance  in  each  waveguide.  This  means  that  the  signal  power  is  decoupled 
from  time  variation  in  the  filter  coefficients. 


0.  Conclusions 

We  have  derived  a  general  framework  for  recursive  digital  filtering  which  has 
many  desirable  properties  (stated  in  the  introduction).  This  architecture  for  digital 
filtering  is  mo6t  valuable  in  the  case  of  time-varying  recursive  networks  in  which  it 
is  desired  to  eliminate  limit  cycles  and  overflow  oscillations.  The  added  complexity 
relative  to  the  best  pre-existing  recursive  filter  architectures  is  negligible.  Therefore, 
these  structures  are  likely  to  find  plentiful  use  in  advanced  time-varying  filtering 
applications. 


7.  Appendix — The  Wave  Equation 

For  an  ideal  waveguide,  we  have  the  following  wave  equation. 

Ptt(z,  t)  -  c2P„(x,  t) ,  (56) 

where  P[z,  t)  denotes  longitudinal  pressure  displacement  in  the  tube  at  the  point  r 
along  the  tube  at  time  t  in  seconds.  If  the  length  of  the  tube  is  L,  then  z  is  taken 
to  lie  between  0  and  L.  The  partial  derivative  notation  used  above  is  defined  by 


The  constant  e  is  given  by  c  =  \fTJ~p  where  T  is  the  tube  “tension,"  and  p  is  the 
mass  per  unit  length  of  the  tube.  An  elegant  derivation  of  the  wave  equation  is 
given  by  Morse  [4j. 

The  general  traveling-wave  solution  to  (56)  is  given  by 

P(z,l)  =  P*(z-et)  +  P~[z  +  ct).  (58) 

This  solution  form  is  interpreted  as  the  sum  of  two  fixed  wave-shapes  traveling  in 
opposite  directions  along  the  tube.  The  specific  waveshapes  are  determined  by  the 
initial  pressure  P(z,  0)  and  flow  V(z,  0)  throughout  the  tube. 
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Abstract 

An  adaptive  equalizer  is  proposed  for  eliminating  distortion  due  to  multipath  propaga* 
tion  in  high-speed  digital  radio  systems.  The  equalizer  is  aimed  at  the  ease  of  very  fast 
channel  fluctuation,  where  ‘‘fast”  is  defined  relative  to  the  impulse-response  duration 
of  the  inverse  of  the  instantaneous  multipath  transfer  function.  The  method  is  of  the 
decision-feedback  type  where  the  demodulated  symbols  are  used  to  construct  an  estimate 
of  multipath-induced  intersymbol  interference.  The  reconstructed  baseband  waveform  is 
then  delayed  and  weighted  according  to  the  current  multipath  parameters,  and  this  simu¬ 
lated  echo  is  subtracted  from  the  incoming  baseband  waveform.  The  distinguishing  feature 
of  our  approach  is  that  an  explicit  model  of  multipath  parameters  replaces  the  tranvsersal 
equalizer  studied  previously. 


1.  Introduction 

Microwave  systems  using  digital  modulation  have  developed  rapidly  since  their 
introduction  in  the  early  1970’s.  This  rapid  growth  is  due  in  part  to  the  tendency 
toward  digital  encoding  of  all  types  of  data,  better  transmission  quality  in  digital 
formats  (especially  in  an  interference  environment),  ease  of  digital  switching,  and 
the  rapidly  falling  costs  of  digital  electronics. 

Multipath  propagation  is  often  the  primary  source  of  error  in  high-speed  digi¬ 
tal  radio,  just  as  it  is  in  the  more  familiar  FM  modulation  systems.  However,  the 
effect  of  multipath  on  microwave-frequency  communication  differs  from  its  effect  on 
analog  FM,  and  different  compensation  techniques  naturally  arise.  For  narrowband 
FM  radio,  outage  is  a  function  of  thermal  noise  or  Sat-fade  margin  [12,2].  For 
digital  transmission,  frequency-selective  fading  causes  severe  amplitude  and  delay 
distortion  which  yield  errors  in  excess  of  Sat-fade  margin  predictions  [5,6].  As  in¬ 
formation  rates  can  be  on  the  order  of  100MHz,  even  small  levels  of  delay  distortion 
can  cause  intersymbol  interference.  For  example,  a  6ns  multipath  delay  can  cause 


complete  digital  eye  closure  in  a  22Mbps  4-QAM  system  operating  at  6GHz,  even 
though  the  received  signal  power  falls  only  lOdB  [12]. 

Techniques  for  minimizing  multipath  distortion  include  frequency  diversity, 
spatial  diversity,  and  equalization.  The  first  two  involve  transmission  on  multiple 
frequencies  or  use  of  multiple  antennas;  when  multipath  degrades  one  frequency 
band  or  antenna  signal,  the  receiver  switches  or  fades  to  another  frequency/antenna. 
Spatial  diversity  alone  can  reduce  multipath  outage  by  a  factor  of  6  to  12  [4]. 
Amplitude  equalization  alone  gives  approximately  a  factor  of  6  improvement  [4]. 
When  spatial  diversity  techniques  are  combined  with  amplitude  equalizers,  the  im¬ 
provement  is  surprisingly  a  factor  of  100  to  800,  because  maximum  power  combiners 
convert  in-band  fade  notches  to  slope,  and  slope  is  what  is  corrected  by  many 
amplitude  equalizers  [12,7].  For  a  minimum-phase  multipath  fade,  slope  equalizers 
often  substantially  correct  the  delay  distortion  as  well  as  the  amplitude  fading. 
However,  for  nonminimum-phase  fades,  a  slope  equalizer  attempts  to  double  the 
phase-delay  rather  than  reduce  it. 

Adaptive  transversal  equalizers  have  been  employed  recently  to  equalize  both 
phase  and  amplitude  [0].  Typically  five  taps  are  used,  where  each  tap  multiplies 
the  received  signal  over  the  time  of  one  symbol.  Thus,  a  linear  combination  of 
five  symbol  periods  is  adaptively  optimized.  The  two  error  criteria  currently  in 
use  are  the  “zero  forcing”  (ZF)  and  “least  mean  square”  (LMS)  errors  [12],  The 
LMS  technique  is  generally  regarded  as  superior  in  performance,  but  more  complex 
to  implement.  An  important  advantage  of  the  transversal  equalizer  is  that  it  is 
equally  effective  against  nonminimum-phase  as  well  as  minimum-phase  fades.  In 
one  study  [10],  outage  reduction  using  spatial  diversity  with  transversal  equalization 
was  a  factor  of  3.3  better  than  that  achieved  using  the  same  spatial  diversity  with 
amplitude  (notch)  equalization  [12]. 

A  further  improvement  can  be  obtained  using  a  so-called  “decision-feedback” 
equalizer  (DFE)  [1,18,14,15,20,21].  Unlike  the  transversal  equalizer,  which  ap¬ 
proximately  inverts  the  multipath  channel  transfer  function,  the  DFE  uses  demodu¬ 
lated  symbols  to  subtract  out  intersymbol  interference  (ISI)  on  later  symbols.  With 
added  processing  delay,  ISI  can  be  subtracted  also  from  earlier  symbols.  An  LMS 
version  of  the  DFE  is  described  in  [1].  The  DFE  is  theoretically  superior  to  linear 
equalizers  in  the  case  of  deep  in-band  notches  [17,18,10].  The  DFE  approach  is 
similar  to  echo  cancellation  [3]  wherein  the  echoes  to  be  canceled  are  constructed 
from  the  estimated  symbol-stream  and  channel  model. 

An  exact  inverse  filter  for  the  multipath  transfer  function  is  recursive,  and  for 
minimum-phase  multipath  an  adaptive  inverse  filter  can  be  utilized  [11,13].  Since 
the  inverse  filter  must  be  constrained  to  a  stable  set,  nonminimum-phase  multipath 
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transfer  functions  cannot  be  inverted  by  a  recursive  equalizer.  However,  under 
certain  conditions  the  multipath  delay  can  be  accurately  estimated  by  tracking  with 
a  stabilized  inverse  filter  even  in  the  non-minimum-phase  case  [13]. 

In  practical  high-speed  digital  radio  systems  (line  of  sight  microwave  links),  we 
have  the  following  characteristics  [12]: 

•  Only  two  or  three  paths  usually  matter. 

•  Delays  range  up  to  approximately  11ns  (multi-GHz  radio). 

•  Distortion  is  due  more  to  intersymbol  interference  caused  by  multipath  delay 
than  by  the  fading  associated  with  multipath. 

•  Fade-notch  depth  — 20  log  1 1  —  |?t[[  can  vary  as  fast  as  lOOdB  per  second, 
where  gt  is  the  multipath  secondary-path  gain. 

«  Notch  frequencies  k/rtt,  k  =  1,2,...,  can  change  as  fast  as  50MHz  per 
second,  where  rt{  is  the  multipath  delay. 

•  Nonminimum-phase  fades  occur  30  to  40  percent  of  the  time,  i .e.,jg^  >  1. 

Under  these  conditions,  it  is  not  very  practical  to  attempt  inversion  of  the 
multipath  transfer  function  H(d,  t). 

This  paper  outlines  a  multipath  equalizer  which  is  based  on  ‘‘decision-directed" 
simulation  of  the  multipath  channel  baseband  output  using  the  demodulated  data 
stream.  The  simulated  channel  output  contains  estimated  echoes  due  to  multipath 
which  are  subtracted  from  the  incoming  baseband  signal.  Our  approach  differs  from 
the  DFE’s  currently  proposed  in  that  the  multipath  is  modeled  explicitly  rather  than 
using  a  transversal  equalizer  to  approximately  model  the  effect  of  multipath  on  the 
baseband.  The  effect  of  hetrodyning  the  multipath  channel  output  is  absorbed  into 
a  complex  multipath  gain  which  is  is  a  function  of  the  multipath  parameters  and 
the  frequency  shift.  All  known  features  of  the  modulation  format  are  thus  exploited 
to  constrain  the  online  optimization. 


2.  Formulation 

Assume  that  the  modulation  format  of  the  transmitted  signal  is  4-QAM  at 
radian  frequency  uc  =  2x/c.  Let  u(t)  denote  the  instantaneous  carrier  frequency, 
Wm  =  2 sr/m  the  modulation  rate,  x(t)  the  baseband  waveform  (the  convolution  of  a 
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rate  um  4- QAM  impulse  train  with  a  Nyquist  shaping  pulse),  and  u{t)  s*  x(f)e/w«< 
the  transmitted  waveform. 

A  model  for  a  time-varying  multipath  channel  is 

H(d,t)=  i +  *</"•  (1) 

where  gt  is  the  gain  of  the  secondary  path  at  time  t,  nt  is  the  instantaneous 
multipath  delay,  and  d  is  the  unit-sample  delay  operator. 

The  received  signal  is  then  v(t)  *  H{d,  t)u(t)  sss  u(t)+gtu(t—nt).  The  recovered 
baseband  signal  is 

y(t)  as  as  *(/)  +  -  ne)  4  x(t)  +  a rtx(f  -  nc)  (2) 

Thus,  the  baseband  signal  appears  as  a  multipath-distorted  signal  itself  where  the 
gain  of  the  secondary  path  is  now  complex: 

(3) 

This  observation  allows  application  of  multipath  cancellation  techniques  which  are 
intended  for  the  modulated  carrier.  (The  modulated  signal  being  in  the  GHz  cannot 
easily  be  dealt  with  directly.) 

The  corresponding  time-varying  inverse  filter  appears  as 

Ht)  **  iKO  -  -  »t) 

»  {z(t)  +  atz{t  -  nt)}  -at{Ht-ht)}  (4) 

*  *(l)  +  o*{x(t  -  n*)  -  x(<  -  ni)} 

The  above  inverse  multipath  filter  is  robust  if  |d<|  <  1  and  the  time  variation  of 
the  multipath  gain  arc  is  slow  compared  with  rtt/(l  —  |aj|).  For  reasons  mentioned 
in  the  introduction,  it  is  not  very  effective  to  attempt  an  exact  linear  inverse  of  the 
multipath  channel. 


3.  Multipath  Equalization 

This  section  describes  a  decision-directed  FIR  equalizer  which  nonlinearly  es¬ 
timates  delayed  copies  of  the  baseband  signal  x(f).  These  multipath  echoes  are  then 
subtracted  out  of  the  received  signal  y(t)  to  give  the  equalizer  output.  A  block: 
diagram  of  the  processing  steps  is  shown  in  Fig.  1. 
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The  input  to  the  system  is  the  received  signal  y(l).  Any  number  of  echoes  may 
be  estimated,  but  for  clarity  we  will  consider  only  one.  Thus,  y{t)  =  x(l)  +  at x(t  — 
fit).  The  estimated  echo  6t[rp(t)\x{t  -  nt[rp(f)J;  r,(<)}  is  subtracted  to  give 

*(0  =  y  (<)  -  d|(rp(<)]x{<  -  ni(rp(/)l;  r.(l)}  (5) 

where  r«(f)  is  the  latest  time  on  which  the  nonlinear  estimate  of  x(t)  is  based,  and 
Tp(t)  b  the  latest  time  used  in  the  parameter  estimates  orj.nj.  The  equalized  signal 
i(t)  is  now  sampled  synchronously  to  recover  the  estimated  information  sequence. 
The  estimated  amplitude  b  quantized  to  the  nearest  signalling  amplitude.  Clearly, 
acceptable  results  will  be  obtained  only  when  thb  quantization  step  produces  the 
true  signalling  amplitude  with  high  probability. 

Given  the  information  sequence,  a  simulated  baseband  it  b  constructed  using 
knowledge  of  the  Nyqubt  shaping  pulse.  Thb  nonlinearly  enhanced  estimate  pos- 
seses  all  known  characteristics  of  the  modulation  format.  The  parameters  of  the 
modulation  format  are  estimated  to  the  extent  that  they  are  unknown;  for  example, 
the  carrier  amplitude  A  and  the  precbe  phase  of  the  switching  transient  would 
normally  need  to  be  estimated  online.  The  basic  idea  b  to  minimize  error  in  the  k- 
step  prediction  of  the  baseband  signal  with  respect  to  switching-time  phase,  carrier 
amplitude,  and  perhaps  even  the  information  sequence.  The  number  k  of  steps  in 
the  prediction  b  determined  by  the  amount  of  processing  delay  in  the  instantaneous 
carrier  estimator,  and  r,(l)  =  t  —  k.  The  mean  square  of  the  difference  it  —  i{t)  b 
minimized  with  respect  to  these  parameters. 

The  multipath  parameters  are  estimated  as  follows.  The  simulated  transmis¬ 
sion  it  b  passed  through  a  simulated  channel 

A.MOI)  - 1 + A.MOI-f*'1''''11  (6) 

to  produce  a  synthetic  version  of  the  received  signal 

y,  *  &(d,(rp(t)|,  n.(rp(0j)x,  (7) 

The  channel  parameters  a t  and  At  are  adaptively  adjusted  to  follow  the  multipath 
gain  at  and  delay  it|.  They  also  must  be  predicted  ahead  of  the  processing  delay.  An 
attractive  criterion  for  minimization  here  b  the  difference  between  the  synthetic 
signal  y<  and  the  received  signal  y(t). 


4.  Description  of  Simulation 


•  A  demodulator  turns  the  equalizer  output  into  a  symbol  stream 

•  The  symbol  stream  drives  a  simulated  transmitter 

•  The  synthesized  transmitter  signal  is  fed  to  the  latest  channel  model 

•  The  echoes  from  the  channel  model  output  are  subtracted  from  the  incoming 
received  signal  to  produce  the  equalizer  output. 

Note  that  the  subtraction  time,  demodulation  time,  synthesis  time,  and  channel- 
model  time  must  all  add  up  to  less  than  the  multipath  delay.  Thus,  the  throughput 
of  the  whole  chain  should  be  on  the  order  of  a  nanosecond,  and  must  therefore  be 
implemented  in  analog  form  (e.g.  optical)  in  VLSI.  The  signal-format  estimation  and 
channel  modeling  proceed  in  parallel  at  slower  rates.  The  signal  format  presumably 
is  close  to  constant.  The  channel  model  changes  at  the  rate  atmospheric  changes 
take  place,  which  is  sufficiently  slow  that  digital  computations  can  be  considered. 

Clearly,  if  the  signal  comes  in  recordable  bursts,  the  delay  from  received  signal 
to  cancellor  signal  is  no  longer  critical,  and  the  multipath  cancellation  can  be  carried 
out  offline. 
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Fig.  1  —  Block  diagram  of  the  proposed  equalizer. 
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TIME-VARYING  AUTOREGRESSIVE  MODELING  OF  A  CLASS 
OF  NONSTATIONARY  SIGNALS 


1.  INTRODUCTION 

The  need  to  estimate  and  process  nonstationary  signals  arises  In  many 
applications.  The  nonstationary  nature  of  the  signals  may  be  caused  by  the 
motion  of  the  signal  source  or  the  receiver  {e.g..  In  radar  or  sonar 
problems),  by  the  properties  of  the  medium  in  which  the  signal  propagates,  or 
by  other  physical  phenomena.  Often  the  signals  are  Inherently  nonstationary 
such  as  In  the  case  of  a  modulated  carrier  In  a  voice  communication  system. 

A  commonly  encountered  type  of  nonstationary  process  consists  of  a 
narrowband  or  sinusoidal  signal  with  a  time-varying  center  frequency  and 
(possibly  stationary)  measurement  noise.  Examples  Include:  the  radar  return 
from  an  accelerating  target,  a  frequency  modulated  carrier,  and  a  variety  of 
acoustic  signals. 

Much  of  the  work  on  the  estimation  and  modeling  of  noisy  signals  Is  based 

on  the  assumption  that  the  signals  may  be  considered  to  be  stationary  over  the 

observation  Interval.  Relatively  little  work  has  been  done  to  take  explicitly 
Into  account  the  effects  of  nonstatlonarlty.  A  promising  approach  to  this 
problem,  which  was  recently  explored  by  several  authors.  Is  to  use  time- 
varying  parametric  models  to  represent  nonstationary  signals.  Linear  models 
of  the  autoregressive  (AR )  and  autoregressive  moving  average  (ARMA)  types  were 
considered  In  [l]-[4].  The  parameters  of  these  models  are  assumed  to  be  time- 

varying,  and  the  functional  dependence  on  time  Is  assumed  to  be  known  up  to  a 

finite  (preferably  small)  number  of  parameters. 

In  this  paper  we  apply  the  time-varying  parametric  modeling  approach  to 
the  class  of  nonstationary  signals  discussed  above.  The  parametric  model  used 
Is  described  in  section  2.  It  Is  shown  that  the  roots  of  the  time-varying  AR 
polynomial,  evaluated  at  a  particular  time  point,  contain  Information  about 
the  Instantaneous  frequency  of  the  signal.  In  a  manner  similar  to  that  of  the 
stationary  case.  In  section  3  we  briefly  summarize  an  algorithm  for 
estimating  the  time-varying  AR  parameters,  which  is  similar  to  the  modified 
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Yule-Walker  (MYW)  method  [5].  Some  examples  Illustrating  the  behavior  of  the 
algorithm  are  presented  in  section  4. 

In  the  stationary  case  the  MYW  equations  provide  an  efficient  (both  in 
the  statistical  and  the  computational  sense)  technique  for  estimating  the  AR 
parameters  of  narrowband  processes.  The  structural  similarity  of  the  MYW 
equations  In  the  stationary  and  nonstationary  cases  may  lead  one  to  believe 
that  the  accuracy  aspects  of  the  resulting  estimates  will  also  be  similar. 
Unfortunately  this  Is  not  the  case.  In  the  nonstationary  case,  estimation 
accuracy  Is  degraded  by  the  presence  of  noise  much  more  than  In  the  stationary 
case.  The  reasons  for  this  behavior  are  briefly  discussed  In  section  3. 
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2.  THE  TIME -VARYING  AUTOREGRESSIVE  MODEL 


We  will  say  that  a  zero-mean  process  x(t)  Is  a  time-varying 
autoregressive  (TVAR)  process  of  order  p.  If  x(t)  obey  the  recursion 

x{t)  «  -  J  a. (t-l)x(t-l)  +  e(t)  (1) 

(So  * 

where  e(t)  Is  a  stationary  white  noise  process  with  zero-mean  and  variance 
2 

a  .  The  time-varying  parameters  {a^ {t),1*l,...,p}  are  assumed  to  be  linear 
combinations  of  a  set  of  basis  time  functions  {f^ft)  ,k«0,. . .  ,m;-  , 

a,(t)  *  I  •  (2) 

i  k«0  IK  K 

The  TVAR  model  Is,  therefore,  completely  specified  by  a  set  of  constant 

2 

parameters  {aik,  1  <  1  <  p,  0  <  k  <  m  ;  0  }  . 

The  choice  of  the  basis  functions  fk(t)  Is  an  Important  part  of  the 
modeling  process.  A  convenient  choice  which  will  be  made  here  Is 

fk(t)  -  (f)k  ,  (3) 

where  T  Is  a  normalizing  constant,  equal  to  the  number  of  data  points  In  the 
observation  Interval.  This  basis  set  Is  by  no  means  the  best  one.  Further 
work  Is  needed  to  develop  a  systematic  procedure  for  designing  optimal  basis 
functions  for  specific  classes  of  nonstationary  signals. 

By  analogy  to  constant  parameter  AR  models  we  can  associate  a  spectral 
function  with  the  autoregressive  parameters.  Let  us  define 

Sx(«;t)  *  *2/|A{z;t)|2,  z-eJ“  ,  (4) 

where 

A(z;t)  ^  l+a1(t-l)z'1+...+ap(t-p)z"p  .  (5) 

The  roots  of  the  polynomial  A{z;t)  can  be  evaluted  for  any  given  time  point. 


The  trajectories  of  these  roots  provide  useful  Information  about  the 
characteristics  of  the  signal  being  modeled. 

To  gain  some  Insight  Into  the  properties  of  the  TVAR  model  consider  the 
simple  case  of  a  sinusoidal  signal  with  linearly  time  varying  frequency, 

x(t)  «  sin  2n(f0+at)t  ,  (6) 

where  fQ  Is  the  Initial  frequency,  and  a  Is  the  frequency  rate  of  change. 

The  Instantaneous  frequency  of  this  signal  Is  given  by, 

f1 (t)  -  fQ  +  2at  .  (7) 

Using  simple  trigonometric  Identities  It  can  be  shown  that  x(t)  obeys 
precisely  the  following  recursion, 

x ( t )  +  b^t-DxU-l)  +  b2(t-2)x(t-2)  •  0,  (8a) 

where 

sin  4»(Mt)+a) 

bl(t*U  *  '  sin  •  ,8b) 

sin  2»(Mt)-a) 

V1-21  *  WTTTT^m^T  •  (8c) 

As  the  parameter  a  tends  to  zero,  these  coefficients  tend 
to  bj  ■  -2  cos  2*fg,  b2  »  1,  which  are  the  parameters  of  the  AR  representation 
of  a  constant  frequency  sinusoid.  More  generally,  for  sufficiently  small 
values  of  a  we  have 

bx (t-1 )  -  2cos  2»Cf0+a(t-l)],  b2 (t-2)  -  1.  (9) 

Thus,  the  roots  of  the  polynomial  A(z;t)  will  be  (approximately)  on  the  unit 
circle  at  angles  ±2irCfg+a(t-l)]  ,  corresponding  to  the  Instantaneous 
frequency  at  time  (t-1).  Evaluation  of  the  spectral  function  Sx(u>;t)  will, 
therefore,  produce  sharp  peaks  at  the  Instantaneous  sinusoidal  frequencies. 


This  analysis  can  be  generalized  to  multiple  sinusoids  as  Is  discussed  In 
[6].  Sinusoidal  signals  whose  frequencies  change  as  a  nonlinear  function  of 
time  can  also  be  approximately  represented  by  TVAR  models.  This  Is  also  true 
for  signals  which  are  not  sinusoidal  but  have  a  non-zero  Instantaneous 
bandwidth.  As  an  example  consider  the  "narrowband"  AR  process  with  time- 
varying  center  frequency  defined  by 

x(t)-2rcos2*f j(t-l)x(t-l)+r2x{t-2)«e(t)  (10) 

where  r  Is  very  close  to  unity. 


3.  THE  MODIFIED  YULE-WALKER  ESTIMATOR 


The  practical  problem  considered  here  Is  one  In  which  It  Is  desired  to 
fit  a  TVAR  model  to  a  finite  number  of  noisy  observations  of  a  signal  of  the 
type  discussed  above.  Let  {y(t),t«l,...,T}be  the  observed  process, 

y(t)  ■  x(t)  +  n(t)  ,  (11) 


where  n(t)  Is  a  noise  process  uncorrelated  with  the  signal  and  x(t)  is  a  TVAR 
process  of  order  p  (1).  It  Is  straightforward  to  show  that 


y(t) 


•*i!i 


a1(t-1)y(t-1)  ♦  v(t) 


(12) 


where 


v(t)  ■  e(t)  +  n(t)  +  y  a.(t-1)n(t-1)  . 

i-i  1 

(13) 

Assuming  that  a^(t)  are  given  by  (2),  we  can  rewrite  (12)  as 

y(t)  ■  -  ?  z( t— 1 )  ♦  v(t) 

1«1  1 

(14a) 

where 

z’(t)  4  (f0(t),...,fm(t))y(t)  , 

(14b) 

A  » 

A1  *  ^a10’***,a1m^  * 

(14c) 

or  In  matrix  notation  (In  the  so-called  autocorrelation  form), 


Let  Z  be  a  shifted  down  and  extended  version  of  the  matrix  Z, 


Z  « 


0 

6 

z*(l) 

• 

z'(T) 

0 


|  q  rows 


N(nH-l)  Columns 


(16) 


The  estimate  of  the  TVAR  parameter  vector  A  Is  computed  by  solving  the  set  of 
overdeterml ned  MYW  equations 

[Z'Z]A  «  -Z'Y  ,  (17) 

where  VI  Is  a  block-toeplltz  matrix  with  N-p  blocks  of  size  (nH-l)x(nH-l) . 

The  estimate  a  can  be  written  as 

A  «  -[Z'ZZ'Z^Z'ZZ’Y  ,  (18) 

although  In  practice  we  would  solve  equation  (17)  without  explicit  matrix 
Inversion. 


To  evalute  the  performance  of  this  estimator  It  is  necessary  to  study  the 
statistics  of  the  estimation  error.  It  Is  straightforward  to  show  that 

A-A  ■  [Z'ZZ'Zj^Z'ZZ’V  .  (estimation  error).  (19) 

In  the  stationary  case  (with  nwl)  It  can  be  shown  that  If  v(t)  Is  finitely 
correlated  and  If  the  delay  parameter  q  (cf.  (16))  Is  appropriately  chosen, 

A 

then  the  error  A-A  will  tend  to  zero  as  the  number  of  data  points  T  tends  to 

A 

Infinity.  In  other  words,  A  In  (18)  Is  a  consistent  estimator  [7].  In  [1] 

It  Is  stated  that  the  estimator  will  also  be  consistent  for  certain  time- 
varying  ARMA  processes. 

The  asymptotic  error  covariance  matrix  for  the  MW  estimates  was 
presented  In  [8], [9]  for  the  stationary  case.  Ho  results  of  this  type  seem  to 
be  available  for  the  nonstationary  case.  However,  some  preliminary  work 
Indicates  that  In  the  case  of  sinusoidal  signals  with  time-varying 
frequencies,  the  estimation  error  of  the  MW  method  may  be  large  compared  to 
the  best  achievable  accuracy  predicted  by  the  Cramer-Rao  bound.  A  heuristic 
explanation  of  this  fact  Is  as  follows:  for  the  MW  method  to  work  well  In 
the  presence  of  noise  It  Is  necessary  that  the  correlation  function  of  the 
signal  decay  relatively  slowly  (cf.  [8]).  In  the  nonstationary  case  the 
entries  of  the  matrix  [Z'Z]  can  be  Interpreted  as  averaged  correlation 
coefficients.  It  can  be  shown  that  for  sinusoidal  signals  with  time  varying 
frequencies  this  "correlation  function"  decays  very  quickly,  leading  to 
Inefficient  parameter  estimates.  In  other  words,  the  MW  estimator  behaves  as 
If  the  signal  was  wideband,  whereas  the  optimal  estimator  makes  full  use  of 
Its  underlying  narrowband  structure.  See  [6]  for  more  details. 


4.  SOME  EXAMPLES 


To  Illustrate  the  behavior  of  the  TVAR  parameter  estimator  we  present  a 
few  examples.  In  all  of  these  examples  the  polynomial  basis  function  In  (3) 
was  used. 

Example  1 

The  signal  consists  of  three  sinusoids  with  linearly  time- varying 
f requencl es 


x(t)*s1n  2»(0.4-8xl0“4t)t+s1n  2»(Q.15+10‘3t)t  (20) 

♦  sin  2w(0.1+8xl0*4t)t 

The  model  order  was  p-6,  and  the  polynomial  order  nw3  (total  of  24 
parameters).  The  number  of  data  points  was  T»128  and  the  signal  was 
practically  noise  free  (SNR»50  dB).  Figure  1  depicts  the  trajectory  of  the 
angles  of  the  roots  of  A(z,t).  Triangles  depict  true  values  and  circles  are 
the  TVAR  estimates.  Note  the  good  fit  achieved  by  the  model,  except  in  the 
neighborhood  of  the  frequency  cross-over  points. 

Example  2: 


Here  the  signal  was  a  single  sinusoid  with  linearly  time-varying 
frequency 


x(t) 


sin  2tt(0.15  + 


I0'4t)t 


(21) 


The  model  order  was  p«2,  the  polynomial  order  was  uwl  and  T-512.  White  noise 
was  added  to  the  signal  to  give  a  signal -to- noise  ratio  of  5  dB.  The  MYW 
equations  were  used  with  q  ■  N  ■  4.  Figure  2  depicts  the  result.  A  good 
model  fit  Is  observed  throughout  the  observation  Interval. 
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Example  3: 

The  TVAR  model  can  be  used  for  signals  with  non- 11  nearly  varying 
frequencies.  In  this  example  we  have  a  single  sinusoid  with  sinusoidally 
varying  frequency, 

x{ t)  »  sin  2*(0.25+0.02  s1n(0.01t))t  .  (22) 

The  model  order  was  p»2,  polynomial  order  nw3,  T*512  and  SNR  >  0  <£.  The  MYW 
equations  were  used  with  q-N>8,  and  the  result  Is  depicted  In  figure  3.  At 
this  low  SNR  the  model  fit  Is  not  very  good  In  some  portions  of  the 
observation  Interval. 

From  these  and  many  other  examples  we  can  make  the  following 
observations:  (1)  the  choice  of  a  model  order  p  equal  to  or  somewhat  larger 
than  twice  the  number  of  sinusoids  yielded  the  best  results.  (11)  the  choice 
of  polynomial  order  m  depends  on  the  maximum  rate  of  change  of  the 
frequencies.  (Ill)  At  low  slgnal-to-nolse  ratios  the  estimator  performed 
poorly,  compared  to  the  optimal  processor. 
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5.  CONCLUSIONS 


Parametric  models  with  time-varying  coefficients  provide  a  powerful 
approach  to  the  representation  and  processing  of  nonstationary  signals.  In 
particular,  the  TVAR  model  was  shown  to  be  useful  for  representing  a  class  of 
narrowband  signals  with  time- varying  center  frequencies.  The  nonstationary 
equivalent  of  the  modified  Yule-Walker  method  was  used  to  estimate  the  TVAR 
parameters.  Some  problems  related  to  the  accuracy  of  this  estimator  were 
briefly  discussed.  It  Is  suggested  that  In  the  nonstationary  case  the 
performance  of  this  estimator  may  be  different  than  expected  from  Its  behavior 
In  the  stationary  case. 
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Abstract 

The  Constant-Modulus  Algorithm  (CMA)  computes  and  applies  an  adaptive 
channel  equalizer  for  constant-amplitude  signals  such  as  frequency-  and  phase- 
modulation.  The  report  “Analysis  and  Extensions  of  the  Constant  Modulus  Algorithm” ' 
by  the  authors  describes  several  extensions  to  the  CMA  Those  which  are  examined 
via  simulations  here  include:  (1)  Extension  from  FIR  equalizers  to  IIR  equalizers 
(poles  as  well  as  zeros  allowed  in  the  equalizer),  and  (2)  Newton’s  method  replaces 
gradient  descent. 


This  work  was  supported  by  the  U.S.  Air  Force  (AFSC),  Route  Air 
Development  Center,  under  Contract  No.  F30602-84-C-0016. 


1.  Introduction 


This  is  an  attachment  to  the  report  “Analysis  and  Extensions  of  the  Constant 
Modulus  Algorithm”  by  the  authors  (hereafter  referred  to  as  the  main  report).  The 
purpose  here  is  to  provide  simulation  results  for  the  algorithms  studied  there.  The 
results  here  are  highly  preliminary. 

For  convenience,  let  CMA-LC  denote  the  LMS  version  of  the  complex  CMA 
algorithm  (gradient  descent),  and  let  CMA-RC  denote  the  RLS  version  (Newton 
descent).  The  simplified  real-only  versions  are  denoted  CMA-LR  and  CMA-RR, 
respectively. 


The  loss  function  minimised  by  the  algorithm  is 


■W-jy  £«?(*) 


where 


.  -  _  i  j|  *t{9)  |  -  mlt,  Complex  Signals 

2 1  x?( 9)  -  0*.,  Simplified  Real  Case 


(1) 

(2) 


1.1.  Gradient  Descent  versus  Newton  Descent 

The  versions  CMA-LC  and  CMA-LR  use  gradient  descent  (essentially  the  LMS 
algorithm),  and  CMA-RC  and  CMA-RR  use  the  recursive  Gauss-Newton  method 
(RGN).  Prior  treatments  of  the  CMA  have  been  based  exclusively  on  gradient 
descent. 

All  four  algorithms  are  implemented  in  their  time- recursive  forms,  that  is,  they 

*  m 

attempt  to  minimize  J(0)  with  respect  the  9  recursively  in  time  as  T  increases. 


The  LMS  version  is 


9t  =  9t-x-liAh-i) 


(3) 


where  is  the  gradient  estimate  of  J{6)  with  respect  to  $  at  time  t,  evaluated 

atf*  0»_i,  and  p  is  a  step-size  parameter  which  is  set  by  the  user. 


The  RGN  version  is  given  by 

(4) 

where  J*(9)  is  the  Hessian,  or  second-derivative  matrix  of  J{6)  with  respect  to  0. 

Before  proceeding  to  the  examples,  we  summarize  the  signal  model  and  algo¬ 
rithm  from  the  main  report. 

1.2.  The  Signal  Model 


Let  yt  denote  the  received  signal  for  1  =  1,2, . . . ,  7*.  In  the  complex  case,  yt  is 
assumed  to  be  of  the  form 


yt  *  +  H%y(d)u t  +  H*i(d)vt  (5) 


where 


(6) 


is  the  transmitted  signal  (y?t  is  the  real-valued  information-bearing  signal),  vt  is 
additive  noise,  ut  is  an  interference  signal  (assumed  at  least  partially  known),  and 
HM9(d)  and  HU9(d)  are  the  linear  time-invariant  channel  filters  associated  with  zt 
and  «i  respectively: 


A(d)  A  do  +  aid  +  fljd3  +  •  •  •  +  antdn* 


B(d)  A  bid  +  63d8  +  ‘'*  +  6n>dB* 
C[d)  A  1  +  Cjd  4*  Cjd2  +  •  •  •  +  Cfi-d1** 


(7) 


where  {a,-,  c,}  are  real  The  unit-delay  operator  d  is  defined  by  the  relation 

A  -«  (8) 


for  an  arbitrary  signal  z%.  We  assume  the  modulus  mx  of  the  transmitted  signal  Z( 
is  known  for  all  t. 

In  the  case  of  real  signals,  the  transmitted  signal  is  assumed  to  be  of  the  form 

zt  *  m,  cos(pt)  (0) 

and  ut  and  vt  are  real  interference  and  additive  noise,  respectively.  We  also  assume 
+  1>t  where 


1.3.  Algorithm  Summary 


Below,  the  notation  c[«]  denotes  the  tth  element  of  the  vector  c. 

y{  =  yt-  - ci-iln*]y/_n4 

(fit  38  (jft,  Sft—i  i ...  i  i  — i  • .  > ,  —  U/— fij  i  —7/-i » *  •  • »  ”1/ 

<p/  ~  I*/,  y/-i . ,  -«Lx . -«£_*, ,  -if.! 

*i  =  pftfi-i 
x!t=*<p{ 

[/?e{x,x,,}1  Complex  Signals 
*i*i,  Simplified  Real  Case 

|2 


s',  *  p(*i*t) 


-  m*,,  Complex  Signals 


x?(0)  —  <rj,,  Simplified  Real  Case 


\  * 


ft,  LMS  version  (Gradient  descent) 
RML  version  (Newton  descent) 


wt 


8tA 


«t(0j,  d([l], . . . ,  d|(n*J,  6,(1], . . . ,  6,[n^j,  c,(l], . . . ,  c,[n*] 
C,(z)-Proj{C,(z)} 

t,  »  tfet 

«**-£,[ l]*^, - c,(nt]Tf_, 

xi{  *  n*  -  c,(l]tt/L1 - ci{n*]iif_, 

*{  *  i,  -  ctfljx/.j - cejntlz/., 


-n» 

•»l 

•"» 


93 


2.  Series  1  —  First-Order  FIR  Equaliser 


In  this  series  of  examples,  the  channel  is  described  by 


H^d) 


1 

So  +  Ojd 


and  the  channel  input  is  given  by 


for  t  a*  0, 1, 2, . . . ,  T  —  1. 


Z|  * 


(12) 


(13) 


The  CMA  estimates  ao  and  aj.  As  shown  in  the  main  report,  in  this  case  the 
CMA  is  guaranteed  to  converge  with  probability  one  to  the  true  solution  [ao,  ai] 
as  long  as  at  least  two  parameters  are  estimated  and  the  phase  modulation  <pt  is 
sufficiently  rich. 


2.1.  Example  1-1  —  Sinsusoidal  FM,  No  Noise 


The  specific  values  used  in  this  example  are 

<pt  am  Uct  +  T0  sin  <Jm 
T  *  2048 
oo  =  1 
ax  —  —0.99 

mx  am  0.1 

P  am  2.375 

u/)n  *  2x0.10 
Ue  am  2x0.11 

94 


a 


Figure  la  shows  the  channel  inpat  signal  Xt,  and  Fig.  lb  shows  the  channel 
output  signal  y<.  The  amplitude  modulation  due  to  the  channel  filtering  is  quite 
visible. 

Figure  2a  shows  the  angle  modulation  waveform,  and  Figure  2b  shows  the 
post-detection  version,  obtained  by  demodulating  y<  =  Hxy(d)xt  to  obtain  <pt. 
Figure  3  gives  the  spectra  of  the  curves  in  Fig.  2.  The  demodulated  signals  have 
some  artificial  discontinuities  due  to  wrap-around. 

Figure  4a  shows  a  close-up  of  the  inverted  normalized  loss  function 

-m~  -i  E  [I *<  Is -mi]’  (U) 

*— 1 

where 

**  **  *©!ft  +  (15) 

and  ao  =  1.  The  peak  is  zero  at  dj  — O.W,  indicating  a  lack  of  bias  in  the 

solution.  Note,  however,  that  the  vertical  axis  coven  only  a  «wiaJ1  range  (near 
the  machine  accuracy  limit);  thus,  the  error  surface  is  quite  flat  near  the  solution, 
indicating  a  large  asymptotic  variance  is  expected  in  the  estimate  of  dj.  This  point 
is  clearer  in  Fig.  4b  which  displays  the  inverted  loss  function  over  a  wider  range  of 

•l- 

Figure  5a  shows  a  close-up  of  the  inverted  normalized  loss  function  (14)  where 
this  time  d\  =  —0.00  and  do  is  varied.  Ti*e  peak  is  zero  at  do  *  1  as  it  must  be. 
Again,  however,  the  error  surface  is  very  flat.  Figure  5b  shows  a  view  from  “farther 
back.* 

A  very  flat  loss  function  means  that  both  the  gradient  J1  and  Hessian  Jn 
are  close  to  zero.  This  means  trouble  for  RGN  which  attempts  to  move  in  the 
direction  J1 .  Newton's  method  is  well-known  to  be  unstable  where  the 

curvature  vanishes.  Root  finders,  for  example,  typically  test  for  this  condition  and 
take  countermeasures  when  necessary.  In  our  case,  a  reasonable  modification  is  to 


add  a  small  constant  times  the  identity  matrix  to  the  Hessian  estimate.  This  is 
accomplished  at  initialization  of  the  algorithm. 

Since  the  third  derivative  vanishes  for  an  even  function,  and  since  a  fourth-order 
Taylor  series  exactly  represents  the  loss  function  (14),  a  descent  method  based  on 
the  fourth-derivative  matrix  (or  tensor)  would  be  expected  to  converge  very  rapidly. 

Figure  6  is  a  3D  plot  of  the  inverted  error  surface  (14)  versus  do  and  ax.  Figure 
7  shows  the  same  plot  from  a  different  angle.  Note  how  flat  the  surface  is  near  the 
middle.  There  almost  appears  to  be  a  ring  of  equal  cost.  However,  we  know  from 
analysis  that  there  are  only  two  local  maxima  in  this  function. 

TheLMSCMA 

First  we  will  run  the  standard  LMS  version  of  the  CMA  to  establish  a  base  for 
comparison.  We  found  empirically  that  ft  =■*  100  was  a  good  setting.  Smaller  ft 
gave  slower  convergence  (by  decreasing  the  slope  magnitude  of  the  essentially  linear 
dx  trajectory).  Larger  ft  was  avoided  because  the  estimate  is  already  starting  to 
oscillate. 

The  initial  conditions  for  this  and  all  later  examples  are  do  *  1,  dx  *  0.  All 
other  initial  state,  when  required  is  initialized  to  zero  unless  otherwise  stated. 

Figure  8  shows  the  estimate  of  dx  produced  by  CMA-LC  for  this  example 
Figure  9a  shows  only  the  first  300  samples  of  the  parameter  trajectory  CMA-LC 
so  that  the  details  of  the  received  signal  yt  (Fig.  9b)  and  the  modulus  error  [|zt|3  — 
m\\f 2  (Fig.  9c)  can  be  seen  in  relationship  to  the  dx  trajectory.  The  final  value  of 
dx  at  sample  2048  is  dx  *  —0.9892  which  is  very  close  to  the  correct  value  =  —0.99. 

Figure  10  shows  the  estimate  of  dx  for  the  real-signal  algorithm  CMA-LR,  and 
Figure  11  gives  the  close-up  view.  We  see  that  the  asymptotic  convergence  is  very 

slow.  In  this  case,  there  is  no  theoretical  assurance  that  the  true  solution  is  obtained 
asymptotically. 


The  RGN  CMA 


Figure  12  shows  the  estimate  of  di  produced  by  CMA-RC  for  this  example. 
Figure  13a  shows  only  the  first  300  samples  of  the  parameter  trajectory  for  CMA- 
RC  so  that  the  details  of  the  received  signal  yt  (Fig.  13b)  and  the  modulus  error 
[|xt|3  —  m3]/ 2  (Fig.  13c)  can  be  seen  in  relationship  to  the  &i  trajectory.  The 
parameter  converges  almost  immediately  to  a  constant  (at  =  —0.9020),  and  the 
modulus  error  rapidly  approaches  a  constant  near  0.175.  While  d|  =  —0.0020  is 
close  to  *  —0.00,  it  is  not  as  close  as  we  expect.  Since  we  know  there  can  be  no 
bias,  it  is  clear  that  the  CMA-RC  algorithm  is  incorrectly  implemented.  That  is, 
we  have  a  “bug*  in  the  simulations  program.  Unfortunately,  time  was  not  available 
to  track  it  down.  Aside  from  the  bias,  note  the  extremely  rapid  convergence  as 
compared  with  the  LMS  version. 

Figure  14  shows  the  estimate  of  a\  for  the  real-signal  algorithm  CMA-RR,  and 
Figure  15  gives  the  close-up  view. 


2.2.  Example  1-2  —  Noise  FM,  No  Noise 


This  example  is  the  same  as  example  1.1  in  every  respect  except  for  the  modula¬ 
tion  signal  which  is  now 

<pt  ■>  Ht(d)t  ,  (16) 

where  c*  is  unit- variance  white  Gaussian  noise,  and  Ht(d)  is  a  6th  order  elliptic 
function  lowpass  filter  with  cut-off  at  one-twentieth  the  sampling  rate  (a  tenth-band 
filter).  The  remaining  figures  are  exactly  analogous  to  example  1.1. 

3.  Conclusions 

The  LMS  version  of  the  CMA  seems  to  be  working  in  both  the  complex  and 
the  real-signal  cases.  The  RML  version  of  the  CMA  seems  to  be  possibly  correct 


■s.  *.  .  ” 


the  real-signal  case,  but  definitely  incorrect  in  the  complex  case. 
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Figure  lb  Spectra  of  Fiaure  1  Curves 
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Figure  5 
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Figure  8 
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Figure  10 
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Figure  11 


''v’V^VW.V  .  V;  ,\V.  V. 

* V  A-  v.j  V%V'V ■  V '  v 


\  «. 


cnun 1C 13 


Tim*  <513:  1  •44-512  (Kf.  Sia:l«34-6l3  ***Of ) 


Clran  1C23  Ti«#  <5  13  :  1*24-5  12  sacs,  512:1*34-513  HWf) 


OSK : *2E . SNO 
OSK : «2H . SNO 


Figure  16 

a)  Angle  modulation  waveform 

b)  Demodulation  of  channel -distorted  signal 
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Figure  19 
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Figure  26 
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Figure  29 
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Figure  31 


i  1 

S  MISSION  i 

5  of  5 

^  Rome  Air  Development  Center  ^ 

R  RAOC  plant  and  executes  Aete.aA.ch,  development,  test  and  v 

*  4  elected  acquisition  program  in  tupport  of  Command,  Control  Q 

6  Communication!,  and  Intelligence  (C3 1)  activities .  Technical  ? 
%  and  engineering  tupport  within  areas  of  technical  competence  C 

it  provided  to  BSV  Program  Offices  { POt )  and  otheA  BSD  > 

5  elements.  The  pAincipal  technical  mittion  aaeat  are  £ 

?  commanicationt,  electromagnetic  guidance  and  control,  tuA-  <* 
^  veillance  o f  ground  and  aero t pace  object. ,  intelligence  data  \ 
V  collection  and  handling,  information  tyttem  technology,  > 

tt  ionotpheric  propagation,  tolid  ztate  tciencet,  microwave 

\  phytic ‘  and  electronic  reliability,  maintainability  and 

6  compatibility. 


